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ABSTRACT 


\ 

The  primary  purpose  of  this  thesis  is  to  to  demonstrate  some  principles  of 
combat  modeling  using  programs  for  the  Radio  Shack  TRS-80  Model  100  computer. 
In  addition  to  the  combat  modeling,  the  thesis  includes  several  utility  programs  for  the 
M100  of  interest  to  students  of  operations  analysis. 

The  combat  modeling  programs  include  an  antisubmarine  warfare  (ASW) 
detection  simulation,  a  Kalman  filter,  and  a  Lanchester  differential  equation 
simulation.  The  utility  programs  include  a  matrix  algebra  program,  a  numerical  double 
integration  program  using  Simpson's  Rule  and  the  Romberg  integration  algorithm,  and 
a  geometric  programming  program  for  zero  degree  of  difficulty  problems.  The 
integration  program  is  also  written  as  a  subroutine  that  can  be  included  in  other 
programs.  The  matrix  algebra  program  includes  a  simultaneous  linear  equation  solving 
subroutine  which  can  be  used  in  other  programs. 

All  programs  are  written  in  M100  BASIC.  Documentation  includes  an 
explanation  of  the  input  required,  the  output  produced,  and  the  components  of  each 
program,  and  sample  problems.  The  chapter  on  geometric  programming  includes  a 
tutorial  on  the  mathematical  basis  for  that  technique. 
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I.  INTRODUCTION 

A.  BACKGROUND 

The  NPS  Operations  Research  Department  has  made  Radio  Shack  TRS-80 
Model  lOO  computers  available  to  a  limited  number  of  students.  Experience  has  shown 
that  these  students  have  made  relatively  limited  use  of  these  computers.  The  main 
reasons  seem  to  be: 

•  This  computer  uses  the  BASIC  programming  language.  NPS  OA  students  do 
not  receive  instruction  in  this  language  but  are  required  to  learn  FORTRAN  and 
APL.  Although  FORTRAN  and  BASIC  have  similar  structures,  most  OA 
students  believe  they  do  not  have  time  to  learn  a  third  computer  language. 

•  Only  a  few  M-100  programs  are  currently  available  at  NPS  that  directly  relate  to 
course  work. 

•  Writing  and  debugging  programs  for  the  M100  can  take  a  considerable  amount 
of  time.  Many  OA  students  believe  that  time  would  be  more  usefullv  spent 
pursuing  other  approaches  to  their  studies. 

•  The  Ml 00  has  thus  far  not  been  distributed  to  all  OA  students.  Therefore, 
professors  have  not  been  able  to  require  students  to  use  the  MlOO.  It  has  been 
relegated  to  "nice  to  have"  status  instead  of  being  included  as  an  essential 
teaching  tool. 

B.  GENERAL 

The  purpose  of  this  paper  is  to  develop  a  collection  of  programs  in  BASIC  for 
the  M 100  that  students  can  use  in  the  combat  modeling  courses  of  the  OA  curriculum. 
The  programs  make  extensive  use  of  subroutines  which  allow  students  to  run  programs 
"off  the  shelf'  or  build  their  own  programs  from  the  subroutines.  The  programs  arc 
extensively  documented  so  that  students  who  learn  BASIC  can  use  the  printed 
programs  as  study  aids  to  understand  the  algorithms  involved.  Some  of  the  programs 
are  tutorial. 

In  testing  situations  professors  are  often  limited  to  problems  that  are 
arithmetically  very  simple.  If  programs  for  the  MlOO  were  available,  professors  would 
be  able  to  give  more  intricate  test  problems  without  placing  undue  emphasis  on  manual 
arithmetic  calculation. 

Additionally,  when  students  leave  NPS,  they  will  be  able  to  take  with  them  a 
series  of  programs  with  which  they  have  grown  familiar  during  their  course  of  study. 
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II.  FEATURES  COMMON  TO  MORE  THAN  ONE  PROGRAM 


A.  GENERAL  CHARACTERISTICS  OF  THE  MODEL  100 

The  M100  is  a  versatile  and  very  portable  computer.  As  provided  to  Naval 
Postgraduate  School  students,  it  has  either  24K  or  32K  of  RAM  for  storage  of 
variables  during  program  execution  and  for  storage  of  programs  and  other  files. 
Programs  and  files  remain  active  while  the  computer  is  turned  off.  There  is  an  internal 
300  baud  modem  which  facilitates  transfering  programs  to  other  computers  for  storage. 
The  eight  line  LCD  screen  limits  the  graphics  display  capability  of  the  MIOO. 
However,  output  in  character  form  can  be  written  to  RAM  files  and  reviewed  after 
program  execution  is  complete.  The  version  of  BASIC  used  in  the  M100,  creation  of 
text  files,  use  of  the  modem,  etc.  are  explained  in  Reference  1.  This  thesis  assumes  that 
readers  are  somewhat  familiar  with  Reference  1. 

Although  the  programs  have  statements  which  print  intermediate  and  final  results 
to  the  screen,  the  operator  may  wish  to  check  the  status  of  a  variable  that  is  not 
printed  by  the  program.  To  do  this 

•  Stop  the  program  by  depressing  the  SHIFT  and  BREAK  keys  simultaneously. 

•  Type  a  BASIC  statement  to  print  the  desired  variables  and  hit  the  ENTER 

button. 

•  Start  the  program  again  at  the  place  it  stopped  by  typing  CONT  and  hitting  the 

EN  I ER  button. 

Most  of  the  programs  in  this  thesis  do  not  include  statements  to  end  the 
program.  This  is  because  the  programs  cycle  back  to  the  start  of  the  program  allowing 
the  operator  to  enter  a  new  set  of  parameters  without  restarting  the  program.  To  end 
the  program 

•  Depress  the  SHIFT  key  and  the  BREAK  key  at  the  same  time. 

•  To  rerun  the  program  from  the  start,  depress  the  F4  key. 

•  To  return  to  the  main  menu,  depress  the  F8  key. 

B.  COMMON  TERMINOLOGY 

BASIC  variables  are  refered  to  in  the  text  using  capital  letters.  Since  M100 
BASIC  only  differentiates  between  variables  based  on  the  first  two  letters  of  the 
variable  name,  most  BASIC  variables  in  the  following  programs  arc  combinations  of 
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two  capital  letters,  e.g.  BA,  CR,  FT.  Names  of  matrices  are  also  specified  using  capital 
letters.  BASIC  permits  matrices  to  be  dimensioned  using  variables.  If,  for  instance, 
AB  =  3  and  AC  =  4,  then  the  BASIC  statement  DIM  ZZ(AB,AC)  would  dimension  a 
three  by  four  matrix  named  ZZ.  When  capital  letters  are  used  inside  the  parentheses, 
the  size  of  the  matrix  is  being  specified.  When  a  small  letter  is  used  inside  the 
parentheses,  a  particular  element  of  the  matrix  is  being  specified.  For  example,  ZZ(i,j) 
refers  to  the  element  of  ZZ  in  the  ith  row  and  jth  column.  When  a  matrix  has  more 
than  one  dimension,  the  first  number  from  the  right  is  refered  to  as  the  column  number 
and  the  second  number  from  the  right  as  the  row  number. 

When  the  mathematical  theory  behind  a  program  is  discussed,  the  variables  used 
will  be  small  letters  or  Greek  letters. 

C.  INPUT  FILES 

All  of  the  programs  presented  in  this  thesis  permit  data  to  be  entered  from  a  text 
file.  These  text  files  must  be  created  using  the  M  100's  text  editor  before  the  program  is 
run.  The  name  of  the  text  file  for  each  program  is  similar  to  the  name  of  the  program 
it  supports  and  is  given  in  the  documentation  for  each  program  under  the  section 
titled,  "Input".  The  contents  of  the  input  file  are  also  explained  in  the  appropriate 
"Input"  section. 

Numbers  in  input  files  must  be  separated  by  a  space,  comma,  or  return.  A 
comma  and  a  return  may  not  be  placed  together  without  a  number  between  them. 
Otherwise^  the  M100  will  enter  an  extra  zero  at  that  point.  A  comma  and  a  return 
together  also  cause  data  elements  which  follow  to  be  in  the  wrong  places  in  their 
matrices  and/or  wrong  numbers  to  be  read  as  matrix  dimensions.  Do  not  put  blank 
lines1  in  the  data  for  the  same  reason. 

D.  FORMULA  TOKENIZATION 

There  are  some  programs  which  must  be  adapted  to  use  different  equations  at 
the  same  point  in  the  program  depending  on  the  application.  For  example,  the 
detection  theory  simulation  in  Chapter  3  must  be  able  to  handle  many  functions  for 
the  probability  of  detection  of  a  sensor.  When  the  function  needs  to  be  changed,  the 
operator  may: 


Two  returns  together. 


iw 


m 


•  Stop  the  program,  call  the  BASIC  editor  for  the  lines  that  need  to  be  changed, 
ana  restart  the  program  after  the  function  has  been  edited,  or 

•  Use  a  tokenization  routine  to  change  the  function  without  stopping  the  program. 
This  section  describes  a  subroutine  (see  Figure  2.1)  which  performs  that  tokenization. 


Figure  2.1  Formula  Tokenization  Section. 

That  subroutine  is  taken  from  [Ref.  2:pp.  58-60].  Tokenization  converts  a  string 
variable  into  an  executable  BASIC  statement.  In  the  example  in  Figure  2.1,  the  right 
hand  side  of  the  function  equation  was  stored  as  a  string  variable,  DFS  before  this 
subroutine  was  called.  For  example,  if  the  equation  was  DF=X  +  Y,  then  DFS  is  the 
string,  "X  +  Y'\  Line  1410  adds  an  equal  sign  and  appropriate  left  hand  side  to  DFS  to 
form  BS,  a  complete  assignment  statement  in  string  format.  In  the  example  above  line 
1410  adds  "DF="  to  DFS  to  form  BS.  Line  1451  computes  memory  location,  Bl,  of 
BS  and  calls  the  machine  level  subroutine  in  the  MI00  read  only  memory  (ROM) 
beginning  at  memory  location  1606.  Subroutine  1606  converts  the  BASIC  keywords  in 
BS  into  their  one  byte  tokens  and  then  stores  them  in  executable  form  at  memory 
location  63105.  Line  1455  calls  the  machine  level  subroutine  beginning  at  memory 
location  2499  which  executes  the  statement  stored  at  memory  location  63105. 


III.  DETECTION  SIMULATION  PROGRAM 


A.  GENERAL 

The  problem  addressed  by  this  program  is  estimating  the  probability  of  detecting 
a  target  whose  location  is  given  by  a  bivariate  normal  distribution  when  the  location  of 
each  detecting  sensor  also  has  a  separate  BVN  distribution.  The  target  distribution  is 
BVN(0,0,<t^  pY  v).  The  location  of  each  sensor,  S;,  is  distributed  BVN(Un  ,  pv  , 

a  y  x,y  1  ul  vi 

<r  ,  <rzv ,  pn  v).  All  distribution  parameters  are  stored  in  an  input  text  file.  The 
ui  i  ui’vi 

program  models  two  sensor  types: 

•  A  Deterministic  Sensor.  This  sensor  has  detection  probabilities  of  1,  0  and  1  in 
three  concentric  circular  bands  around  the  sensor's  location.  An  example  of  this 
type  of  sensor  would  be  a  sonobouy  which  detects  all  targets  out  to  range  r  j , 
detects  no  targets  between  ranges  rj  and  detects  all  targets  in  a  convergence 
zone  between  Tj  an^  r-j,  and  detects  no  targets  beyond  r-j  where  0  <  rj  ^  ^  ^ 

r3- 

•  A  Probabilistic  Sensor.  This  sensor  has  a  continuous,  radially  symmetric 
detection  pattern.  The  probability  of  detection,  Pd,  is  a  function  of  range  and 
may  also  be  a  function  of  one  or  more  other  parameters.  The  default  function, 
D(r,b),  for  calculating  Pd  is  a  Carlton  function  where  b  is  a  scaling  parameter 
(See  Equation  3.1).  Figure  3.1  shows  graphs  of  several  Carlton  detection 
functions. 

2  2 

Pd  =  D(r,b)  =  e'r  !2b  (eqn  3.1) 


For  either  type  of  sensor  the  program  calculates  the  overall  probability  of  detection 
using  D(r).  The  deterministic  sensor  portion  of  the  program  can  also  be  used  to 
calculate  a  "cookie  cutter"  approximation  of  the  actual  detection  function.  The 
cookie-cutter  approximation  has  the  form  D(r)=  1  when  r^rQ  and  D(r)  =  0  when 
r>rQ,  for  some  specified  detection  range,  rQ.  The  radius,  r,  of  the  cookie-cutter 
approximation  is  the  lethal  radius  of  the  actual  detection  function  for  a  sensor.  Lethal 
radius  is  a  term  borrowed  from  the  damage  functions  of  firing  theory.  It  is  used  here 
to  help  readers  who  are  familiar  with  firing  theory,  not  because  being  detected  by  a 
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The  program  uses  a  Monte  Carlo  simulation  to  estimate  P^.  Because  the  Monte 
Carlo  simulation  is  time  consuming,  faster  numerical  integration  approximations  are 
also  available  for  the  deterministic  sensor.  However,  the  particular  numerical  integration 
technique  used  in  this  program  will  not  produce  correct  results  if  sensor  areas  of  coverage 
overlap.  More  sophisticated  numerical  integration  techniques  are  available  for  special 
cases  of  the  probabilistic  sensor,  but  are  not  included  in  this  program. 

B.  EXPLANATION  OF  VARIABLES 

•  AL  is  the  probability  that  the  computed  confidence  level  does  not  contain  the 
true  probability  of  detection  for  the  associated  standard  normal  CDF  value. 

•  A2(6,6)  is  the  intermediate  calculation  matrix  for  the  Romberg  integration 
subroutine. 

•  BO  and  B1  are  the  memory  location  parameters  of  BS. 

•  BS  is  the  string  that  holds  the  entire  equation  for  DF  in  the  tokenization  process. 

•  D1  and  D2  are  intermediate  calculation  values. 

•  DF  is  the  value  of  the  detection  function. 

•  DFS  is  the  string  that  hold  the  right  hand  side  of  the  equation  for  DF. 

•  F  is  the  value  of  the  function  being  integrated  at  a  specific  point. 

•  FI  is,  in  the  data  input  section  (lines  100- 150),  the  selection  variable  indicating 
whether  all  sensors  have  the  same  parameters.  If  FI  =  1,  then  sensor  parameters, 
other  than  the  five  for  location,  wiil  be  listed  only  once:  after  the  location 
parameters  of  Sj.  If  FI  =  0,  then  values  of  all  parameters  for  all  sensors  must  be 
entered  explicitly  in  DSIN.DO.  In  the  rest  of  the  program,  beyond  line  200,  FI 
is  a  flag  determining  whether  the  calculation  is  based  on  the  deterministic 
(FI  =  1)  or  the  probabilistic  (FI  =  2)  sensor. 

•  F2  is  the  selection  variable  for  whether  all  sensors  locations  are  distributed  with 
<rx  =  <Ty  =  0.  1  -  yes;  2  -  no. 

•  F3  is  the  selection  variable  for  whether  a  Monte  Carlo  simulation  (F3=  1)  or  a 
numerical  integration  (F3  =  2)  is  used. 

•  H  is  the  radius  of  circular  limits  of  integration. 

•  I,  II,  12  are  loop  counters. 

•  IN  is  the  value  of  the  integral  from  a  numerical  approximation. 

•  J1  through  J9  are  loop  counters. 

•  NS  is  the  number  of  sensors. 

•  NP  is  the  number  of  parameters  for  each  sensor  in  addition  to  location. 


•  NR  is  the  number  of  repetitions  in  a  Monte  Carlo  simulation. 

•  PD  is  the  probability  of  detection. 

•  RFis(l-RH2)'5. 

•  RH  is  the  correlation  between  components  of  the  BVN  target  location 
distribution. 

•  SI  and  S2  are  standard  deviations  of  the  components  of  the  target  location 
distribution. 

•  TE  is  a  temporary  storage  variable. 

•  T1S  and  T2S  are  times  at  beginning  and  end  of  a  calculation. 

•  T1  is  the  time  elapsed  during  a  calculation. 

•  VI  and  V2  are  variances  of  the  components  of  the  target  location  distribution. 

•  XS  and  XT  are  specific  values  of  the  first  component,  i.e.  the  mean  X  position, 
of  the  BVN  distributions  of  a  sensor  or  the  target  respectively. 

•  X(NS,5  +  NP)  is  the  parameter  matrix  for  sensors.  There  is  one  row  per  sensor 
and  one  column  per  parameter.  Columns  one  through  five  are  for  the  means, 
u„  o  and  liv  o  ,  the  standard  deviations,  <T.,  Q  and  c  ,  and  the  correlation, 

U,3j  V,0-  U,Sj  V,5j 

PS:’  of  the  location  of  each  sensor,  S-.  These  parameters  are  in  the  same 
coordinate  system  as  the  target  location  distribution. 

•  YS  and  YT  are  specific  values  of  the  second  component,  i.e.  the  mean  Y  position, 
of  the  BVN  distributions  of  a  sensor  or  the  target  respectively. 

•  Z1  and  Z9  are  selection  variables. 

C.  INPUT 

1.  Input  File 

Before  this  program  is  run,  a  text  file,  DSIN.DO,  must  be  created  to  hold  the 
input  parameters.  DSIN.DO  will  contain  the  following  variables  in  the  following 
order: 

•  NS,  NP,  SI,  S2,  RH,  FI 

•  The  sensor  location/parameter(s)  matrix,  X(NS,N'P  +  5). 

An  example  of  a  input  file  is  shown  in  Figure  3.2  along  with  a  graph  of  the 
corresponding  sensor  coverage  areas.  Entries  in  the  first  line  of  Figure  3.2  indicate  that 
there  are  two  sensors  with  three  parameters  each,  that  the  target  location  distribution 
is  BVN(0,  0,  625,  625,  0),  and  that  all  sensors  have  the  same  parameters.  The  second 
line  is  the  parameter  set  for  the  first  sensor.  There  is  no  variance  in  the  location  of 


Figure  3.2  Example  of  Input  File.  DSIN, 
And  Graph  Of’ Sensor  Coverage. 


either  sensor.  The  first  sensor  is  located  at  (20,20)  and  has  additional  parameters  5,  25, 
and  30.  If  this  input  file  represents  sensors  that  have  a  deterministic  detection 
functions,  then  the  sensors  detect  all  targets  within  a  distance  of  5  and  within  a  circular 
band  from  25  to  30  and  do  not  detect  targets  outside  of  these  bands.  Both  sensors 
have  the  same  parameters  other  than  location,  as  indicated  by  the  last  element  of  the 
first  line  being  one,  the  last  line  contains  only  the  location  of  the  second  sensor.  If  one 
or  more  sensors  had  different  parameters  (other  than  location),  then  the  last  element  in 
the  first  row  would  have  been  zero.  Also,  the  three  nonlocation  parameters  would 
have  been  specified  for  all  sensors.  Although  the  parameters  for  each  sensor  may  be 
different,  the  detection  function  for  each  sensor  must  be  the  same.  For  example,  the 
program  does  not  allow  a  mixture  of  sensors  with  deterministic  and  probabilistic 
detection  functions. 

2.  Interactive  Input 

After  the  program  has  read  the  input  file,  the  operator  will  be  prompted  to 

specify: 

•  Whether  a  deterministic  or  a  probabilistic  detection  function  will  be  used. 

•  Whether  the  variance  in  target  location  is  or  is  not  equal  to  zero. 

•  Whether  a  Monte  Carlo  simulation  or  numerical  approximation  will  be  used. 
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If  the  probabilistic  detection  function  is  used,  the  operator  may  interactively  change 
the  detection  function.  If  the  Monte  Carlo  simulation  is  used,  the  operator  will  be 
prompted  to  specify  the  number  of  repetitions,  NR,  in  the  simulation. 

D.  OUTPUT 

The  program  prints  the  probability  of  detection,  PD,  to  the  screen  of  the  M100. 
If  the  computation  was  a  Monte  Carlo  simulation,  the  program  asks  the  operator  to 
specify  an  a3  for  the  confidence  interval  on  the  probability  of  detection.  Permissible 
a  s  are  .10,  .05,  and  .01.  The  program  computes  the  lower  and  upper  limits  of  the 
confidence  interval  based  on  a  normal  approximation  to  the  binomial  distribution.  The 
expression  for  these  limits  is  shown  in  Equation  3.3. 

PD  ±  Z\.a  CPD(1-PD)/NR]-5  (eqn  3.3) 

/2 

Zj_a  is  the  point  on  a  standard  normal  distribution  at  which  the  CDF  has  a  value  of 
l -a  2  2 

E.  EXPLANATION  OF  PROGRAM  COMPONENTS 

A  complete  program  listing  is  at  Appendix  A. 

1.  Initialization/Data  Input,  Figure  3.3 

Lines  100-115  print  the  program  title  and  open  the  input  file,  DSIN. 

Line  120  reads  the  first  line  of  DSIN  and  calculates  the  variances  of  the 
marginal  distributions  of  the  target  location  distribution.  It  also  dimensions 
X(NS,5  +  NP)  and  A2(6,6). 

Lines  125-135  read  the  lines  of  DSIN  that  specify  parameters  for  individual 
sensors,  S- ,  i=  1,2,. ..,NS.  Line  125  reads  all  parameters  for  Sj.  If  some  Sj  have 
different  nonlocation  parameters,  i.e.  F1=0,  then  line  12S  reads  in  5  +  NP  parameters 
for  Sj  where  i=2,...,NS.  If  all  sensors  have  the  same  nonlocation  parameters,  i.e. 
FI  =  1,  then  lines  132-135  read  the  five  location  parameters  for  each  sensor.  They  also 
make  nonlocation  parameters  of  Sj,  where  i  =  2,3,4,...,NS,  equal  to  the  nonlocation 
parameters  of  S j. 


3a  is  the  probability  that  the  computed  confidence  level  does  not  contain  the  true 
probability  of  detection. 


100  CLS : PRINT" " : PRINT"  DETECTION  SIMULATION" : FORI =1T0400 :NEXTI 

110  'Input/Initialization 
115  0PEN"DSIN“F0RINPUTAS1 

120  INPUT#1 ,NS ,NP ,S1 ,S2 ,RH , FI : V1=S1*S1 : V2=S2*S2 

121  DIMX<NS,5+NP>,A2(6,6),T1<3) 

125  F0RI1=1T05+NP : INPUT#1 ,X<  1 ,11 ) :  NEXTIl 

126  IFNS=1THEN140 

127  IFF1=1THEN132 

128  F0RI1=2T0NS: F0RI2=1T05+NP :  INPUTR1  ,Xt  II ,12  ) :NEXTI2 :NEXTI1 :G0T0140 
132  F0RI1=2T0NS: F0RI2=1T05; INPUT81 »X(  II ,12 ) :NEXTI2:IFNP=0THEN135 

134  F0RI2=6T05+NP:X( 11,12 )=X(  1 ,12  >  :NEXTI2 

135  NEXTIl 

140  RF=SQR<  1-RHa2) 

150  DF$="FXP(  -( ( XT-XS  )A2+(  YT-YS  )A2  )/(  2*X(  J2 ,6  )A2  ) )" 


Figure  3.3  Initialization  and  Data  Input  Section. 


Line  140  calculates  RF,  an  intermediate  value  used  in  later  integration 
calculations. 

Line  150  assigns  the  right  hand  side  of  the  Carlton  equation  to  DFS  as  the 
default  equation  for  calculating  the  detection  function  value,  DF. 

2.  Method  Selection  Section,  Figure  3.4 


200  '*»*  Simulation  Selection  Section  **» 

201  CLS:PRINT"Is  the  Detection  function:” 

203  PRINT"  1.  Deterministic": PRINT"  2.  Probabilistic" 

205  INPUT-Enter  1  or  2:"jFl 

210  CLS:PRINT"Are  Sensor  Locations:" 

212  PRINT"  1.  Always  At  Aim  Point": PRINT"  2.  Distributed  BVN  Around  Aim  Point" 

214  INPUT"Enter  1  or  2:">F2 

215  IF( F1=20RF2=2  )THENF3=1 :G0T0230 

220  CLS: PRINT"" :PRINT"Is  the  Calculation:" 

222  PRINT"  l=Monte  Carlo  Simulation" : PRINT"  2=Numerical  Approximation" 

224  INPUT"Enter  1  or  2"»F3 

230  TIME$="00: 00:00": IF  F1=1THENGOSU8300ELSEGOSUB500 
250  GQT0200 


Figure  3.4  Method  Selection  Section. 

Lines  200-250  assign  values  to: 

•  FI.  If  F 1  =  1 ,  then  calculations  are  done  for  a  deterministic  sensor.  If  FI  =  2, 
then  calculations  are  done  for  a  probabilistic  sensor. 
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•  F2.  If  F2=  1,  then  the  sensor  is  always  located  at  its  aiming  point,  i.e. 

cr^x  =  =  0.  If  F2  =  2,  then  cr^x  or  G^y  is  not  equal  to  0. 

•  F3.  If  F3=l,  then  a  Monte  Carlo  simulation  is  conducted.  If  F3  =  2,  then  a 
numerical  approximation  is  calculated. 

Line  230  sets  the  M  100's  clock  to  zero  to  time  the  calculation.  It  also  branches  to  the 
appropriate  subroutine  depending  upon  which  type  of  sensor  is  specified  by  FI. 

3.  Deterministic  Sensors,  Lines  300-360 


300  'Deterministic  Sensor  Subroutine 

305  IFF3=1THENGOSUB310ELSEGOSUB350 

306  RETURN 


Figure  3.5  Decision  Logic  For  Deterministic  Sensors. 


Based  upon  the  value  of  F3,  lines  300-306  (see  Figure  3.5)  branch  to 
subroutines  to  conduct  the  calculations  using  a  Monte  Carlo  simulation  or  a  numerical 
approximation. 

a.  Actual  Detection  Function ,  Lines  310-360 

(1)  Monte  Carlo  Simulation,  Figure  3.6.  Line  315  calls  subroutine  900 
which  prompts  the  operator  for  the  number  of  repetitions,  NR. 


310  'Monte  Carlo  of  Deterministic  Sensor 
315  GOSUB900 

320  PD-0 : F0RJ1=1T0NR : PRINT . 241  /'Repetition: ” ) J1 : GOSUB600 : FORJ2=lTONS 

323  IFF2=2THENXS=X( J2,1):YS=X<  J2 ,2 ) :GOT0325 

324  GOSUB612 

325  Tl=SQRt ( XS-XT I A  2  +  <  YS-YT ) A  2  ) 

330  IFT1<=X< J2,6  )THENPD=PD+1 ; GOT0335 

332  IFT1>=X<  J2 , 7 !ANDT1<=X( J2 ,8 )THENPD=PD+1 : GOT0335 

334  NEXTJ2 

335  NEXTJ1:PD=PO/NR:GOSUB950: RETURN 


the  circular  detection  bands  around  each  sensor,  S-,  i=  1,2 NS.  PD  counts  the 

number  of  repetitions  for  which  the  target  is  in  at  least  one  detection  band.  The 
Monte  Carlo  simulation  functions  accurately  even  if  detection  bands  of  various  sensors 
overlap  because  it  does  not  double  count  if  a  target  is  detected  by  more  than  one 
sensor.  When  all  repetitions  are  completed,  the  counter,  PD,  is  divided  by  NR  to 
produce  an  estimate  of  the  probability  of  detection.  Finally,  output  subroutine  950  is 
called. 

(2)  Numerical  Approximation,  Figure  3.7.  The  volume4  under  the  target 
location  distribution  is  calculated  for  the  detection  bands  of  each  S^.  This  volume  is 
the  sensor's  probability  of  detection.  Overlapping  detection  bands  are  not  permitted  with 
this  numerical  approximation  because  the  program  would  double  count  the  overlapping 
volumes.  Subroutine  1200  conducts  the  integration.  Probability  of  detection  is  the 
volume  under  the  target  location  distribution  that  is  within  the  coverage  area  of  any  Sj. 


350  ‘Numeric/Oeterminis'tic  Subroutine 

355  PD=0:FORJ2=1TONS:H=X( J2 ,6  ) : GOSUB12O0 : PD=PD+IN 

356  H=X<  J2 ,8  ) : GOSUB1200 : PD=PD+IN 

357  H=X<  J2 , 7  ) : GOSUB1200 : PD=P0-IN:NEXTJ2 
360  GOSUB950 : RETURN 


Figure  3.7  Numerical  Approximation. 

4.  The  Probabilistic  Sensor,  Figure  3.8 

Lines  502-503  print  a  header  and  call  subroutine  1300  which  permits  the 
detection  function  to  be  modified.  Line  503  also  calls  the  subroutine  in  which  the 
operator  specifics  the  number  of  Monte  Carlo  repetitions. 

For  each  repetition,  lines  520-530  generate  a  target  location  from  the  BVN 
target  location  distribution.  Then,  for  each  S-  they  calculate  the  value  of  the  sensor's 
detection  function  using  subroutine  1410  and  generate  a  uniform  random  number 
between  zero  and  one.  If  that  random  number  is  less  than  or  equal  to  Sj's  detection 
function  value  at  that  location,  then  Sj  detected  the  target  and  PD  is  augmented  by 


4  This  volume,  although  it  is  computed  and  described  as  a  volume  in  this  section, 
is  not  a  true  volume.  This  is  because  although  the  X  and  Y  axes  of  the  BVN  are 
distances,  the  Z  axis  is  a  probability,  i.e.  dimensionless.  Therefore  the  result  is  really 
an  area,  not  a  volume. 
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500  'Probabilistic  Detection  Function 

502  CLS:PRINT"De'fault  Detection  Function  Is  Carleton." 

505  GOSUB1300:GOSUB900 

520  P0*0 : F0RJ1=1T0NR : PRINT . 241 , "Repetition:" > J1 : GOSUB600 : F0RJ2=1T0NS 

521  IFF2=2THENXS=X<  J2,l ) : YS=X<  J2,2  )  .-G0T0523 

522  GOSUB612 

523  GOSU81410:IFRNS( 1  )<=DFTHENPD=PD+1 :G0T0526 

524  NEXTJ2 

526  NEXTJ1:PD=PD/NR 
530  GOSU8950: RETURN 


Figure  3.8  The  Probabilistic  Detection  Function. 


one.  If  one  sensor  detects  the  target,  the  program  moves  directly  to  the  next, 
repetition.  Therefore,  the  Monte  Carlo  simulation  functions  accurately  even  if  detection 
bands  of  various  sensors  overlap  because  it  does  not  double  count  a  target  that  is  detected 
by  more  than  one  sensor.  When  all  repetitions  are  completed,  the  counter,  PD,  is 
divided  by  NR  to  produce  an  estimate  of  the  probability  of  detection.  Finally,  output 
subroutine  950  is  called. 

5.  Generating  A  BVN  Random  Variable,  Lines  600-606 


600  1 ***Genera  to  BVN  RV*** 

602  U1=RND( 1 )  :U2=RND( 1 ) :TE=SQRl -2*L0G(U1 ) ) 

604  XT =TE*COS(  6 . 2831853*U2  ) :  YT=RH*XT+RF*TE*SIN<  6 . 2831853*112  ) 

606  XT=XT*S1:YT=YT*S2: RETURN 

612  U1=RNDI 1 )  :U2=RND( 1 )  :TE=SQR( -2*LOG( U1 ) ) 

614  XS=TE*COS<  6 . 2831853*U2 ) : YS=X( J2 ,5 )*XT+(  1-X(  J2 ,5 )*2 ) A . 5*TE*SIN<  6 . 2831853*U2 ) 
616  XS=X<  J2 , 1 )+XS*X( J2 ,3 ) : YS=X<  J2 , 2 )+YS*XI J2 ,4 ) : RETURN 


Figure  3.9  Generating  A  BVN  Random  Variable. 


Lines  602-606  generate  the  target  location.  Lines  602-604  compute  both 
components  of  a  BVN(0,  0,  1,  1,  0)  random  variable  using  two  independent  uniform 
(0,1)  random  numbers  generated  by  the  M  100's  RND(l)  function.  liquations  3.4,  and 
3.5,  as  described  in  [Ref.  4:page  953],  form  the  basis  for  lines  602  and  604.  Equations 
3.4  generate  two  independent  normal(0,l)  random  variables,  Xj'  and  X2'. 


Xj'  =  (-21n  Uj)'5  cos(2rtU2)  (cqn  3.4) 

X2'  =  (-21n  U,)-5  sin(27tU2) 


26 


Equations  3.5  convert  these  two  independent  normal  random  variables  into  the 
components  of  a  BVN(0,  0,  1,  1,  p)  distribution,  Xj  and  X2. 

Xj  =  Xj'  (eqn  3.5) 

X2  =  pXj'  +  (l-p2)-5X2 

Line  606  scales  the  components  of  the  BVN(0,  0,  1,  1,  px  )  by  SI  and  S2  to 
produce  the  components  of  the  BVN(0,  0,  <T  x,  cr  y,  px  y)  target  location  distribution 
and  then  returns  to  the  calling  program. 

Lines  612-616  generate  a  sensor  location  using  the  same  algorithm  that  was 
used  to  generate  the  target  location.  However,  in  addition  to  being  scaled  by  (Tu  and 
Gy,  sensor  locations  must  also  be  displaced  by  pu  and  py. 

6.  Entering  The  Number  Of  Repetitions,  Figure  3.10 


Figure  3.10  Entering  The  Number  Of  Repetitions. 

Line  900  prompts  the  operator  to  interactively  specify  the  number  of 
repetitions,  NR,  for  Monte  Carlo  simulations. 

Line  910  is  a  subroutine  which  stops  the  program  while  the  operator  views 
output  to  the  screen. 

7.  Output  Section,  Figure  3. 1 1 

All  printing  is  to  the  screen  of  the  M100.  Lines  951-952  play  a  tune  to  notify 
the  operator  that  the  calculation  is  finished.  Line  953  also  prints  the  time  required  to 
do  the  calculation  and  branches  around  the  selection  of  a  for  the  confidence  interval  if 
a  numerical  approximation  was  used.  Lines  954-959  prompt  the  operator  to  select  .1, 
.05,  or  .01  as  a  =  AL,  the  probability  that  the  true  probability  of  detection  is  not  in  the 
confidence  interval.  Once  AL  is  selected,  it  is  reassigned  the  value  of  the  standard 
normal  at  Zj_tt/2.  Line  960  prints  the  point  estimate  of  Pj.  Line  962  prints  a  short 
reminder  that  there  is  no  confidence  interval  with  numerical  approximations.  Line  966 
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950  'Print  output 

951  SOUND 1567 > 10 :S01M)1244, 10: SOUND1046  >10 : SOUND 783 >20 

952  SOUND1046,10:SOUND783,40 

953  CLS:  PRINT'"' :  PRINT  "Calculation  Time  (HH/MM/SS)  =  "jTIME*:IFF3=2THEN960 

954  PRlNT"S«lect  Alpha  for  Confidence  interval:'* 

955  INPUT"  Choices  =  .1,  .05,  .01:"»AL 

956  IFAL= . 1THENAL=1 . 645 : GOTO960 

957  IFAL= . 05THENAL=1 . 96 : 60T0960 

958  IFAL= . 01THENAL=2 . 575 : GOT0960 

959  GOT 0954 

960  PRINT'****  Estimate  of  P<Detection)  =  "» :PRINTUSING"  #.#####"»PD 

961  IFF3=1THEN965 

962  PRINT“No  Confidence  Interval  For  Numerical  Approximations" 

963  GOTO970 

965  PRINT"Conf i dance  Interval:  ("> 

966  TE=AL*SQR( P0*( 1-PD  )/NR ) : LL=PD-TE :UL=PD+TE :IFUL>1THENUL=1 

967  IFLL<0THENLL=0 

968  PRINTUSING"## . #####" » LL jUL  j :  PRINT"  )*' :  GOSUB910 
970  'Confetti  Approximation 

972  PRINT"" :INPUT"Confetti  approximation?  0=No,  l=Yes:"»Z9:IF29=0THENRETURN 
974  CLS: INPUT"Enter  TOTAL  lethal  area  for  ALL  sensors  in  the  pattern :"jNA 

976  TE=NA/(6.283185*S1*S2 ):TE=l-( l+SQRI 2*TE ) )*EXP< -SQR( 2*TE ) ) 

977  PRINT“**Confetti  Approximation  =  ")TE :GOSUB910: RETURN 


Figure  3.11  Output  Section. 

calculates  the  confidence  interval  according  to  Equation  3.3.  Lines  966-967  ensure  that 
the  confidence  limits  are  between  zero  and  one.  Lines  965  and  968  print  the  confidence 
interval. 

Lines  970-977  approximate  Pd  with  the  "confetti  approximation"  described  in 
[ReE  5:pp.  14-16].  This  approximates  Pd  by  distributing  the  total  lethal  area  of  the 
entire  group  of  sensors  over  an  ellipse  as  described  in  Reference  5.  This  approximation 
is  calculated  by  Equation  3.6  where  £  =  (Total  Lethal  Area)/(27tcrx<Ty). 


Pd=  1-(1  +  (2Q-5)  eW 


(eqn  3.6) 


Line  972  prompts  the  operator  to  indicate  whether  a  confetti  approximation  is 
desired  and  ends  the  output  routine  if  it  is  not.  Line  974  prompts  the  operator  to  enter 
the  total  lethal  area  for  all  sensors  combined.  Line  976  computes  C,  and  then  Pj.  Line 
977  prints  Pd  and  ends  the  output  subroutine. 

8.  Numerical  Integration,  Figure  3.12 

This  subroutine  is  a  tailored  version  of  the  Romberg  integration  subroutine 
documented  in  Chapter  8.  It  integrates  the  target  location  distribution  subject  to 
circular  limits  of  integration. 


1200  ’Numerical  Integration  Subroutine 

1201  Dl=6. 2831853*S1*S2*RF 

1202  TLS.001 

1220  C  LS : PRINT "": PRINT "  ! 'Calculating  An  Integral!  !,•:PRINT,,■, 

1230  N=2:GOSUB1293:DY=( YU-YL )/2 

1240  FORJ9=1T06:OY=DY/2:N=N*2 

1242  Y=YU: GOSUB 1296: GOSUB1 280 : A2 ( J  9 , 1 )  =TS*DX 

1245  Y=YL : GOSUB1296 : GOSUB1280 : A2( J9 ,1  )=A2<  J9, 1  )+TS»DX 

1250  F0RJ8=2T0N: Y=Y+DY :G0SUB1296 : GOSUB1280 

1251  A2( J9,l  )=A2( J9,l  )+2*TS*DX:NEXTJ8 

1252  A2(J9,1)=A2<  J9,l)*DY/2 
1255  IFJ9=1THENNEXTJ9 

1260  F0RJ8=1T0J9-1 

1262  A2(  J9,J8+1)=A2(  J9,J8)  +  (  <A2<  J9,J8)-A2(  J9-1.J8)  )/(4AJ8-l)  ):NEXTJ8 

1263  T1=A2<  J9,  J9  )-A2<  J9,  J9-1 ) : IFSGNI T1  )*T1-TL>0THENNEXTJ9ELSE1266 

1264  PRINT"Tolerance  of"jTLj"not  met  after  five  extrapolations" 

1266  IN=A2U9,J9): RETURN 

1275  F0RJ7=1T06 : F0RJ6=1T0J7 : PRINTUSING"## . ###" l A2( J7 , J6  ) » 

1276  NEXTJ6 : PRINT"" : NEXTJ7 : INPUTZ9 : RETURN 

1280  REM  Trapezoidal  Rule  Sum 

1281  X=XU : GOSUB 1286 : TS=F : X=XL : GOSUB1286 : TS=TS+F 

1282  F0RJ5=2T0N-1 : X=X+DX : G0SUB1286 : TS=TS+F : NEXT J5 : RETURN 

1285  '  F(X,Y)  to  be  integrated: 

1286  F=Xa2/V1-2*RH*X*Y/S1/S2+Ya2/V2 

1287  F=(EXP(-F/2/RFA  2  n/01:  RETURN 
1290  ‘Limits  of  Integration: 

1293  YU=X(J2,2)+H:YL=X( J2, 2 )-H: RETURN 

1296  T3=SGR(Ha2-(Y-X(  J2,2 ) )a2  ) :XU=X(  J2,l  )+T3 :XL=X(  J2,l )-T3:0X=(XU-XL)/N 

1297  RETURN 


Figure  3.12  Numerical  Integration  Subroutine. 
9.  Changing  The  Detection  Function,  Figure  3.13 


location.  XS  and  YS  represent  the  location  of  Sj2  in  the  same  coordinate  system  as 
XT  and  YT.  An  example  of  a  probabilistic  detection  function  that  is  not  Carlton 
might  be  an  exponential  function  as  specified  in  Equation  3.7. 


Pj  =  Xe‘^r  where  r  =  distance  from  sensor  to  target.  (eqn  3.7) 

The  entry  to  be  made  when  prompted  by  line  1 306  would  be 

X(J2,6)*EXP(-X(J2,6)*SQR((XT-XS)A2  +  (YT-YS)A2)) 

where  X  for  sensor  J2  is  X(J2,6).  In  general,  X(J2,6),...,X(J2,5  +  NP)  represent  NP 
other  parameters  of  the  selected  detection  function  in  addition  to  the  five  parameters  of 
the  BVN  sensor  location  distribution.  The  formula  for  the  current  detection  function  is 
displayed  by  line  1304.  The  operator  may  either  change  it  as  desired  or  hit  ENTER  to 
leave  the  current  formula  unchanged. 

10.  Formula  Tokenization  Section,  Figure  3.14 


1400  'Tokenize  OF 

1410  B$="DF="+0F$+CHR$( 0  ) 

1450  ' T oken iza/exGcu te  B$ 

1451  B0=VARPTR( B$  ) : B1=PEEK( BO+1 )+256*PEEK< B0*2 ) : CALL1606 ,0 ,B1 
1455  CALL2499 , 0  > 63105 : RETURN 


Figure  3.14  Formula  Tokenization  Section. 

The  right  hand  side  of  the  detection  function  equation  is  stored  as  a  string 
variable,  DFS.  This  section  converts  DFS  into  an  executcable  BASIC  assignment 
statement  and  executes  that  statement.  A  detailed  explanation  of  this  subroutine  is 
found  in  the  Formula  Tokenization  Section  in  Chapter  2. 

F.  SAMPLE  PROBLEMS 

Examples  1-6  which  follow  have  the  following  characteristics  in  common. 

•  They  use  a  target  location  distribution  that  is  BVN'(0,  0,  25  ,  25  ,  0),  except  for 
Example  6  in  which  pY  =  .5. 

•  All  measurements  are  in  kilometers  and  are  in  the  coordinate  system  of  the  target 
location  distribution. 

•  All  confidence  intervals  are  calculated  with  a  =  .05. 

30 


•  All  Monte  Carlo  simulations  were  done  with  5,000  repetitions. 

•  All  sensors,  deterministic  or  probabilistic,  have  a  lethal  radius  of  twenty  and 
therefore  a  lethal  area  of  40071. 

1.  Example  1 

This  example  is  of  a  single  cookie-cutter  sensor  located  at  (0,0)  with  a  lethal 
radius  of  20.  The  input  file  and  diagram  of  sensor  coverage  is  at  Figure  3.15.  The 
closed  form  solution  may  be  calculated  using  Equation  3.8. 


1  -  e-202;(2*252)  .  .273s3 


(eqn  3.S) 


1  3  25  25  0  0 
00000  20  00 


»  SENSOR  COVERAGE  (SHADED  RRERS) 

X  TGT  LOG  IS  OIST  BVN1D.0, 525,625,0! 
LOCATION  DISTRIBUTION  (DOTTED  LINES) 


Figure  3.15  Input  File  And  Sensor  Coverage  Diagram  For  Example  1. 

The  probabilities  of  detection  computed  by  the  various  computational 
techniques  and  calculation  times  were  as  follows. 

•  Closed  form:  .27385 

•  Numerical  integration:  .27252;  3  minutes,  24  seconds 

•  Monte  Carlo  simulation:  .27640  ±  .0124;  4  hours,  22  minutes. 
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2.  Example  2 

This  example  is  of  a  single  Carlton  sensor  located  at  (0,0)  with  a  lethal  area 
equivalent  to  the  cookie-cutter  sensor  in  Example  1;  i.e.,  b=20/(2’5)  =  14.14214.  The 
input  file  and  diagram  of  sensor  coverage  is  shown  in  [Ref.  5:page  5).  The  closed  form 
solution  is  at  Equation  3.9.  Because  the  sensor  is  located  at  the  origin  and  <r  =  <j 
Equation  3.9  simplifies  to  Equation  3.10  for  this  example. 


Figure  3.16  Input  File  And  Sensor  Coverage  Diagram  For  Example  2. 

The  probabilities  of  detection  computed  by  the  various  computational 
techniques  and  calculation  times  were  as  follows. 

•  Closed  form:  .242424 

•  Numerical  Integration:  None 

•  Monte  Carlo  simulation:  .2452  ±  .01 19;  3  hours,  44  minutes. 


Pj  =  y  e''^u  +  ^v) 

where  y  =  b2/  [(b2  -f  <j2x)(b2  +  <r2v)] 

6u  =  M2u/(b2  +  <*2U).  and  5v  =  M2v.,(b2  +  ff2v)- 


(eqn  3.9) 


(eqn  3.10) 


Pd  =  b2  (b2  +  <r2)  =  14.1422/(14.1422  +  252)  =  .242424 
3.  Example  3 

This  example  is  for  a  single  cookie-cutter  sensor  offset  at  {10,10)  with  a  lethal 
radius  of  20.  The  input  file  and  diagram  of  sensor  coverage  is  shown  in  Figure  3.17. 
There  is  no  closed  form  solution  for  offset  cookie-cutter  sensors. 


1  3  25  25  0  0 

10  10  0  0  0  20  0  0  *  SENSOR  COVERAGE  (SHADED  AREAS) 

■  TGI  LOC  IS  OIST  8VNI0, 0,625.625,0) 
LOCATION  DISTRIBUTION  (DOTTED  LINES) 


Figure  3.17  Input  File  And  Sensor  Coverage  Diagram  For  Example  3. 

The  probabilities  of  detection  computed  by  the  various  computational 
techniques  and  calculation  times  were  as  follows. 

•  Closed  form:  None. 

•  Numerical  integration:  .2400;  2  minutes,  51  seconds. 

•  Monte  Carlo  simulation:  .24160  ±  .01 187;  2  hours,  57  minutes. 

4.  Example  4 

This  example  is  for  a  single  Carlton  scnsoi  offset  at  (10,10)  with  a  lethal  area 
equivalent  to  the  cookie-cutter  sensor  in  Example  3,  i.e.  b=20/(2"^)=  14.14214  Tnc 
input  file  and  diagram  of  sensor  coverage  is  shown  in  Figure  3.18.  The  equation  for 
the  closed  form  solution  is  shown  in  Equation  3.9. 


Figure  3.18  Input  File  And  Sensor  Coverage  Diagram  For  Example  4. 


The  probabilities  of  detection  computed  by  the  various  computational 
techniques  and  calculation  times  were  as  follows. 


•  Closed  form:  .21475 

•  Numerical  integration:  None 

•  Monte  Carlo  simulation:  .2138  ±  .0114;  3  hours,  15  minutes. 

5.  Example  5 

This  example  is  for  a  single  convergence  zone  sensor  iocated  (0,0)  with  a  lethal 
area  equivalent  to  the  sensors  in  Examples  1-4.  The  sensor  detects  all  targets  at  ranges 
less  than  four,  detects  no  targets  between  four  and  30,  detects  all  targets  between  30 
and  35.83,  and  detects  no  targets  beyond  35.83.  The  input  file  and  sensor  coverage 
diagram  is  shown  in  Figure  3.19.  The  closed  form  solution  is  found  by  using  Equation 
3.8  to  calculate  Pd  in  circles  of  radii  4,  30,  and  35.83.  Then  Pd  =  Pd(4)  -  P^f  30)  + 
Pd(35.83). 

The  probabilities  of  detection  computed  by  the  various  computational 
techniques  and  calculation  times  were  as  follows. 

•  Closed  form:  .141463 

•  Numerical  integration:  .14128 

•  Monte  Carlo  simulation:  .1416  ±  .0097;  2  hours,  59  minutes. 
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Figure  3.19  Input  File  And  Sensor  Coverage  Diagram  For  Example  5. 

6.  Example  6 

This  example  is  for  four  convergence  zone  sensors  with  mean  locations  at 
(10,10),  (-10,-10),  ( 10.- 10),  and  (-10,10).  These  sensors  have  convergence  zones  equal 
to  that  of  the  the  sensor  in  Example  5.  Sensor  locations  are  distributed  BVN  with  pu 
and  pv  as  indicated  above  and  <ru=<Tv=3,  and  puy  =  .7.  The  input  file  is  at  Figure 
3.20.  There  is  no  closed  form  solution  for  this  example. 


The  probabilities  of  detection  computed  by  the  various  computational 
techniques  and  calculation  times  were  as  follows. 

•  Closed  form:  None 

•  Numerical  integration:  None 
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IV.  KALMAN  FILTER  PROGRAM 


A.  GENERAL 

This  program  successively  updates  an  estimate  of  the  state  of  a  system,  p,  based 
upon  a  series  of  measurements,  Vectors  p  and  £  need  not  have  the  same  dimensions. 
The  program  makes  provision  for  changes  in  the  system  state  between  measurements  in 
accordance  with  a  linear  system  model.  Covariances  between  elements  of  p  and 
between  elements  of  C,  are  also  accounted  for  by  the  program. 

The  purpose  of  this  chapter  is  to  describe  an  implementation  of  a  Kalman  filter 
on  the  M100.  A  fuller  explanation  of  the  mathematical  background  behind  Kalman 
filters  may  be  found  in  References  6  and  7.  The  notation  in  the  following  program  is 
similar  to  the  notation  used  in  Reference  6  to  facilitate  comparison  of  the 
mathematical  theory  and  computer  implementation. 

The  state  of  the  system,  X,  is  assumed  to  be  a  multivariate  normal  random 
variable,  X  ~  N(p,E),  with  system  noise,  W  ~  (pw,Q),  and  measurement  noise,  V  ~ 
(pv,R).  (p  is  the  matrix  which  models  the  linear  change  in  X  between  measurements. 
H  is  the  matrix  which  shows  the  linear  relationship  between  the  measurement  and  the 
system  state,  i.e.  how  the  measurement  depends  on  the  system  state. 

The  Kalman  filter  recursively  updates  p  by  repeating  two  steps,  measurement  and 
movement.  The  measurement  step  calculates  the  Kalman  gain,  enters  a  new 
measurement,  and  updates  p  and  L  based  on  that  measurement.  The  movement  step 
updates  p  and  L  based  upon  the  system  model. 

B.  EXPLANATION  OF  VARIABLES 

•  B 1  is  the  intermediate  matrix  for  inversion  of  C2. 

•  C1(MD,MD)  and  C2(MD,MD)  are  the  matrices  which  hold  intermediate  results 
of  matrix  calculations. 

•  II  is  the  matrix  showing  the  relationship  between  measurements  and  the  system 
state. 

•  II,  12,  and  13  are  loop  counters. 

•  K(NX,NZ)  is  the  Kalman  gain  matrix. 

•  MD  is  the  maximum  of  NX  and  XZ. 

•  ML’(NX)  is  the  current  estimate  of  the  system  state. 

•  MV(NZ)  is  the  mean  of  measurement  noise. 

•  MW(NX)  is  the  mean  of  system  noise. 
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•  NX  is  the  number  of  elements  in  the  system  state  vector. 

•  NZ  is  the  number  of  elements  in  the  measurement  vector. 


•  PH(NX,NX)  is  the  matrix,  <p,  modeling  the  linear  change  in  the  svstem  state 
between  measurements. 

•  Q(NX,NX)  is  the  covariance  matrix  of  system  noise. 

•  R(NZ,NZ)  is  the  covariance  matrix  of  measurement  noise. 

•  SG(NX,NX)  is  the  covariance  matrix  of  the  system  state. 

•  Z(NZ)  is  the  vector  holding  the  current  measurement,  £. 

C.  INPUT 

The  operator  must  create  a  RAM  file,  KALIN. DO,  which  contains  the  program 
input.  KALIN. DO  must  contain  the  following  parameters  and  matrices  in  order. 

•  NX  and  NZ,  the  number  of  elements  in  ft  and  C,  respectively. 

•  Matrix  PH(NX,NX),  the  system  model  matrix,  (p. 

•  Vector  MW(NX),  the  mean  of  system  noise,  ^w. 

•  Matrix  Q(NX,NX),  the  covariance  matrix  of  system  noise. 

•  Matrix  H(NZ,NX),  the  matrix  showing  the  relationship  between  ft  and 

•  Vector  MV(NZ),  the  mean  of  measurement  noise,  ny. 

•  Matrix  R(NZ,NZ),  the  covariance  matrix  of  measurement  noise. 

•  Vector  ML’(NX),  the  initial  estimate  of  the  system  state. 

•  Matrix  SG(NX,NX),  the  initial  estimate  of  the  covariance  matrix  of  the  svstem 
state,  E. 

After  the  initial  estimate  of  the  system  state  is  entered  from  the  input  file, 
KALIN. DO,  the  program  will  prompt  the  operator  for  measurements  and  changes  to 
the  H  matrix  to  be  entered  from  the  keyboard. 

D.  OUTPUT 

All  output  goes  to  the  screen  of  the  M 100.  After  the  measurement  step  the 
program  displays  the  updated  Kalman  gain  matrix,  the  estimated  system  state,  fi  +  , 
and  the  covariance  matrix,  E  + .  After  the  movement  step  the  program  displays  the 
estimate  of  the  system  state,  ft-,  and  covariance  matrix,  E-,  just  prior  to  the  next 
measurement. 

E.  EXPLANATION  OF  PROGRAM 

A  complete  program  listing  is  located  at  Appendix  B. 


1.  Initialization/Input,  Figure  4.1 


100  CLS:PRINT"***#*KALMAN  FILTER)*#***"  :  PRINT"  Input  Data  Being  Read" 

110  0PEN"KALIN"F0RINPUTAS1 : ONERRORGOT09900 
120  INPUT81,NX,NZ:IFNX<NZTHENMD=NZELSEMD=NX 

125  DIMPHI NX, NX  )  ,MH( NX  I  ,Q( NX, NX  )  ,Hl NZ.NX  )  ,MV( NZ  )  ,Rl NZ.NZ  )  ,MU( NX  )  ,SG( NX, NX  ) 

126  DIMC1(MD,MD  ),C2(MD,MD  ) ,K( NX ,NZ  )  ,B1( NZ+1 ,NZ*2  ) 

130  F0RI1=1T0NX : F0RI2=1T0NX : INPUTS 1 >PH(I1,I2)  :NEXTI2 :NEXTI1 
132  F0RI1=1T0NX:INPUT(?1,MW(I1  l-.NEXTU 

134  F0RI1=1T0NX:F0RI2=1T0NX:INPUT81,Q< II ,12  ) : NEXTI2 : NEXTI1 
136  F0RI1=1T0NZ : F0RI2=1T0NX: INPUT81 , H( II ,12 ) : NEXTI2 : NEXT II 
138  F0RI1  =  1T0NZ : INPUT81 ,MV( II  ):NEXTI1 

140  F0RI1=1T0NZ : F0RI2=1T0NZ: INPUT81 ,RII1,I2):NEXTI2: NEXTI1 
142  F0RI1=1T0NX : INPUTffl ,MU( II): NEXTI1 

144  FORI 1=1T0NX : F0RI2  =  1T0NX : INPUT81 ,SG( II ,12  ) : NEXTI2 : NEXTI1 

145  00=0 : CLS : PRINT " In i tial  SG  As  Input  Check GOSUB532 


Figure  4.1  Initialization  and  Input  Section. 

Line  110  opens  the  input  file,  KALIN. DO  and  branches  the  program  to  line 
9900  if  an  error  occures.  Line  120  enters  NX  and  NZ  and  calculates  MD,  Lines  125 
and  126  dimension  the  matrices  in  the  program.  Lines  130-144  enter  the  matrices  from 
KALIN. DO  as  described  in  the  input  section  above.  Line  145  initializes  the 
measurement  counter,  CC,  and  prints  the  last  input  matrix,  SG,  as  a  check  on  input 
accuracy. 

2.  Measurement  Block,  Lines  150-387 

Measurement  block  matrix  equations  for  updating  K,  ji,  and  L  are  listed  in 
Equations  4.1,  4.2,  4.3. 

K  =  LIlVllLII1  +  R)'1  (eqn  4.1) 


M  =  B  +  K(Z  -  Hv  -  lift)  (eqn  4.2) 


E  =  (I  -  KII)E 


(eqn  4.3) 


Sffl. 

•M 


*#1 


c.  Entering  A  New  H  Matrix ,  Figure  4.2 


ISO  CLS: PRINT"  *****MEASUREMENT  BLOCK*****" 

160  PRINT "Cur rent  H  ":GOSUB540 

162  INPUT"Enter  New  H  ?  l=Yes,  0=No:"»Z9:IFZ9=0THEN170 
165  ‘Enter  A  New  H 

167  F0RI1=1T0NZ:F0RI2=1T0NX 

168  PRINT”Enter  Row"(Ill">  Column" »I2 1 “Of  H 

169  INPUTHI 11,12 ) :NEXTI2 : PRINT"" :NEXTI1 


Figure  4.2  Entering  a  New  H  Matrix. 

Some  Kalman  filter  problems  require  a  different  H  matrix  for  each 
measurement.  This  section  permits  such  a  matrix  to  be  entered.  Line  160  prints  a 
header  and  calls  the  subroutine  which  prints  the  current  H  matrix.  Line  162  asks  the 
operator  whether  a  new  H  matrix  is  required  and  branches  the  program  appropriately. 
Lines  167-169  prompt  the  operator  to  fill  the  new  H  matrix  row  by  row. 
b.  Kalman  Gain  Calculation ,  Figure  4.3 


170  ‘CALC  KALMAN  GAIN 

171  ‘MULT  SG  H  t,  INTO  Cl 

172  FORIl=lTONX: FORI2=1TONZ:C1( 11 ,12  1=0 : F0RI3=1T0NX 

174  Cl( 11,12  )  =  (SG( 11,13  )*H( 12,13  ))+C II 11,12  ):NEXTI3:NEXTI2:NEXTI1 
180  ‘MULT  H  SG  H  1  i  INTO  C2 

182  F0RI1=1T0NZ: F0RI2=1T0NZ:C2( 11,12 )=0: F0RI3=1T0NX 

184  C2(I1,I2)  =  (H(I1,I3  )*C1(  13, 121) +C  2(11, 12) :  NEXTI3 :  NEXTI2 :  NEXTI1 

200  ‘ADO  R  INTO  C2 

202  FORI 1=1T0NZ: FORI 2=1T0NZ:C2( 11,12 )=C2I 11,12 )+R( 11,12) 

203  NEXTI2:NEXTI1 
210  'INVERT  C2 
215  GOSUB9800 

220  ‘MULT  Cl  C2  INTO  K 

222  FORI1=1TONX:FORI2=1TONZ:K(I1,I2)=0:FORI3=1TONZ 

224  K(I1,I2)  =  (C1(I1,I3)*C2(I3,I2))+K(I1,I2):NEXTI3:NEXTI2:NEXTI1 


Figure  4.3  Kalman  Gain  Calculation. 

This  section  updates  the  Kalman  gain  matrix,  K,  in  accordance  with 
Equation  4.1.  Lines  171-174  multiply  E  by  II1  and  place  the  result  in  matrix  Cl. 
Lines  180-184  multiply  H  by  EH1  and  place  the  result  in  matrix  C2.  Lines  200-203  add 
R  to  (HEM1)  and  place  the  result  in  C2.  Line  215  calls  the  subroutine  which  inverts 
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(HEH*  +  R).  Lines  220-224  multiply  EH1  by  (HEH1  +  R)  ,  producing  an  updated 
Kalman  gain  matrix,  K. 

c.  Enter  A  Measurement  and  Update  p,  Figure  4.4 


250  '***##UPDATE  MU-  TO  MU+«*** 

251  'MULT  H  MU-  INTO  Cl 

252  F0RI1=1T0NZ : C1C II ,1 )=0 : F0RI3=1T0NX 

254  C7.<  11,1  )=IH(  11,13  )*MU<  13)  )+Cl< II ,1 ) :NEXTI3:NEXTI1 
260  'ADO  MV  +  H  MU- 

262  F0RI1=1T0NZ:C1< 11,1 )=C1( 11,1 )+MV< II )  :NEXTI1 
270  'INPUT  A  NEW  MEASUREMENT 

272  CC=CC+1 : CLS : PRINT "Measurement  #")CC)":" 

273  FORI 1=1T0NZ : PRINT "Enter  Element" jll»"0f  Measurement:") 

274  INPUTZI II ) :NEXTI1 

280  'SUBTRACT  Cl  FROM  Z,  INTO  Cl 

282  F0RI1=1T0NZ:C1( 11,1 )=Z(  II )-Cl(  11,1 )  :NEXTI1 

290  'MULT  K  Cl  INTO  C2 

292  F0RI1=1T0NX:C2( 11,1 )=0 : F0RI3=1T0NZ 

294  C2<  11,1  )  =  (  Kill,  13  )*C1<  13, 1 )  )+C2l  11,1)  :NEXTI3  :NEXTI1 

300  'ADO  C2  ♦  MU-  TO  UPDATE  TO  MU+ 

302  F0RI1=1T0NX:MU( II )=C21 11,1 )+MU(  II ):NEXTI1 


Figure  4.4  Enter  Measurement  And  Update  The  Estimate  of  ft. 

This  section  enters  a  new  measurement,  Z,  from  the  keyboard  and  updates 
the  estimate  of  in  accordance  with  Equation  4.2.  Lines  251-254  multiply  LI  by  p  and 
place  the  result  in  Cl.  Line  262  adds  py  +  Hp  and  places  the  result  in  Cl.  Lines 
272-274  increment  the  measurement  counter  and  allow  the  operator  to  enter  the  new 
measurement,  Z.  Line  282  subtracts  (py  +  Hp)  from  Z,  and  places  the  result  in  Cl. 
Lines  292-294  multiply  the  Kalman  gain,  K,  by  (Z  -  py  -  Up),  and  place  the  result  in 
C2.  Line  302  adds  p  to  K(Z  -  py  -  Hp),  producing  the  revised  estimate  of  p. 

d.  Updating  E,  Figure  4.5 

This  section  updates  the  estimate  of  E  using  Equation  4.3.  Lines  322-326 
multiply  the  Kalman  gain,  K,  by  H,  subtract  the  result  from  the  identity  matrix,  I,  and 
place  the  result  in  CL  Lines  352-354  multiply  (1  -  KM)  by  E  and  place  the  result  into 
C2.  This  result,  the  updated  E,  is  then  copied  into  SG  by  line  362. 

e.  Printing  The  Updated  K,  p,  And  E,  Figure  4.6 

Lines  375-377  print  a  header  and  call  the  printing  subroutine  for  the 
updated  Kalman  gain.  Lines  380-382  print  a  header  and  call  the  printing  subroutine 
for  the  updated  estimate  of  the  system  state,  ft.  Lines  385-387  print  a  header  and  call 
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320  ’MULT  K  H  S  SUBTR  FROM  I  ,  PUT  IN  Cl 

322  FORI 1=1T0NX: FORI 2=1T0NX: Cl  ( 11,12  )=0 : F0RI3=1T0NZ 

324  Cl(  11,12  )  =  (K(  11,13  )*H(  13,12  )  )+Cl(  II  ,12  ):NEXTI3:C1(  11,12  )  =  -Cl<  11,12) 
326  NEXTI2:NEXTI1 

328  F0RI1=1T0NX:C1( 11,11 )=1+C1< 11,11  ):NEXTI1 

3S0  'MULT  LAST  RESULT  BY  SG  ,  INTO  C2 

352  FORI 1= 1TONX : FORI 2=1T0NX : C2( II ,12 )=0 : F0RI3=1T0NX 

354  C2( 11,12 )=(C1( 11,13  )*SG( 13,12  )  )+C2( II ,12 > :NEXTI3:NEXTI2:NEXTI1 

360  'PUT  C2  INTO  SG 

362  FORI1=1TONX: FORI 2=1T0NX: SGI  II, 12) =C 2(11, 12  )  :NEXTI2:NEXTI1 


Figure  4.5  Updating  The  System  Covariance  Matrix,  E. 


375  CLS:PRINT"Kalman  Gain,  Kli.jJ  After" 

377  PRINT "Measurement  «"jCC:GOSUB510 

380  CLS:PRINT“Estimate  Of  System  State,  MUC i  )+  After" 

382  PRINT "Measurement  #"lCC:GOSUB520 

385  CLS:PRINT"Estimate  Of  Covar,  SG<i,j)  +  After" 

387  PRINT"Measurement  #"  »CC:G0SUB530 


Figure  4.6  Printing  The  Updated  Kalman  Gain,  p,  And  E. 

the  printing  subroutine  for  the  updated  estimate  of  covariance  matrix  of  the  system 
state,  E. 

3.  Movement  Block,  Lines  400-490 

Movement  block  matrix  equations  for  updating  p  and  E  are  listed  in 
Equations  4.4,  and  4.5. 


P  =  <PP  + 


E  =  tpEtp1  +  Q 


(cqn  4.4) 


(eqn  4.5) 


a.  Updating  The  Estimate  of  p,  Figure  4.7 

Lines  422-424  multiply  <p  by  p  and  place  the  result  in  Cl.  Line  432  adds 
pw  to  <pp,  producing  the  updated  estimate  of  p  just  before  measurement  CC+  1. 


a 


a 


$ 

M 


400  CLS:PRINT"*********MOVEMENT  BLOCK*******" 

410  'Update  MU(CC)  +  to  MUICC+1)- 

420  'MULT  PH  MU  ,  PUT  IN  Cl 

422  F0RI1=1T0NX:C1<I1,1 1=0: F0RI3=1T0NX 

424  C1<I1,1)  =  <PHII1,I3)*MU(I3))+C1<I1,1):NEXTI3:NEXTI1 

430  'ADD  C1+  MW  ,  INTO  MU 

432  FORIl=lTONX:MU( II )=C1< 11,1 >+MW( II  ):NEXTI1 


Figure  4.7  Updating  The  Estimate  of  p. 


b.  Updating  E,  Figure  4.8 


440  '**UPDATE  SG  ** 

4S0  'MULT  PH  SG  ,  INTO  Cl 

452  F0RI1=1T0NX : F0RI2=1T0NX:C1(  II ,12 )=0 : F0RI3  =  1T0NX 

454  C1(I1,I2)  =  (PH(I1,I3  )*SG(  13, 12)) +C 1(11, 12  ):  NEXTI 3  :  NEXTI 2 :  NEXTI  1 

460  'MULT  Cl  PH  t,  INTO  C2 

462  F0RI1= 1TONX : FORI 2 = 1TONX :C2(I1,I2)=0:FORI3 = 1TONX 

464  C2(  11,12  )  =  (C1(  11,13  )*PH(  12,13)  )+C  2111, 12  )  :NEXTI3 :  NEXTI  2 :  NEXTI  1 

470  'ADD  C2  ♦  Q  =  SG 

472  F0RI1=1T0NX: FORI2=1TONX:SG< 11,12 )=C2< 11,12 )♦£)(  11,12 ) 

473  NEXTI 2: NEXTI 1 


Figure  4.8  Updating  E. 


Lines  452-454  multiply  (p  by  E  and  place  the  result  in  Cl.  Lines  462-464 
multiply  <pE  by  <p*  and  place  the  result  in  C2.  Lines  472-473  add  Q  to  <pE(p\ 
producing  the  updated  estimate  for  E  just  before  measurement  CC+  1. 

c.  Printing  The  Updated  p  And  E,  Figure  4.9 


480  PRINT"Estimate  Of  System  State,  MU( I 
482  PRINT"Before  Measurement  If"  )CC  +  1:GOSUB520 
485  CLS:PRINT"Estimate  Of  Covar,  SG(I,J)-  Before" 
487  PRINT"Measurement  ff"fCC*l:GOSUB530 
490  GOTO160 


Figure  4.9  Printing  The  Updated  p  And  E. 
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Lines  480-482  print  a  header  and  call  the  printing  subroutine  for  the 
updated  estimate  of  the  system  state,  p.  Lines  485-487  print  a  header  and  call  the 
printing  subroutine  for  the  updated  estimate  of  the  system  state  covariance  matrix,  L. 
Line  490  branches  the  program  back  to  the  beginning  of  the  measurement  block. 

4.  Printing  Subroutines,  Figure  4.10 

Lines  500-554  contain  subroutines  which  print  the  Kalman  gain  matrix,  ft,  £, 
the  H  matrix,  or  the  C2  matrix. 


500  'PRINTING  SUBROUTINES 
510  'PRINT  KALMAN  GAIN,  K 

512  FORIl=lTONX: FORI2=lTONZ: PRINTUSING"######. ##"jK< II, 12 )j:NEXTI2 
514  PRINT'"' : NEXTI 1  :INPUT"H it  ENTER  To  Continue:  "  JZ9: RETURN 
520  'PRINT  MU 

522  FORI 1=1T0NX: PRINTUSING"######.##" vMU( II ) i :NEXTI1 : PRINT"" 

524  INPUT''Hit  ENTER  To  Continue »Z9:RETURN 
530  'PRINT  COVAR  MATRIX,  SG 

532  FORIl=lTONX: FORI 2=1T0NX : PRINTUSING"#####. ##" >SG( 11,12 )> : NEXTI 2 
534  PRINT"" : NEXTI 1 : INPUT"H it  ENTER  To  Continue:"jZ9:RETURN 
540  'PRINT  H 

542  FORI 1=1T0NZ : FORI 2= 1TONX : PRINTUSING"###### . ##" >  H  ( 1 1 , 1 2 ) J : NEXTI 2 
544  PRINT"" :NEXTI1: RETURN 
550  PRINT"  C2  MATRIX:" 

552  FORIl=lTOA: FORI 2=lTOB: PRINTUSING"###### . ##" >C2(  II ,12 ) > : NEXTI 2 
554  PRINT"" :NEXTIl:INPUT"Hit  ENTER  To  Continue:"jZ9:RETURN 


Figure  4.10  Printing  Subroutines. 

5.  Inversion  Subroutine  For  C2,  Figure  4.11 

This  subroutine  is  essentially  the  same  as  the  matrix  inversion  subroutine  in 
the  matrix  algebra  program  discussed  in  Appendix  E.  It  has  been  abbreviated  to  invert 
only  matrix  C2  instead  of  inverting  several  matrices  of  various  dimensions  as  in 
Appendix  E.  This  subroutine  also  does  not  calculate  the  determinant  of  C2  to  test  for 
invertability.  If  C2  is  not  invertablc,  an  division  by  zero  error  will  occur  and  the 
program  will  branch  to  the  error  identification  section.  A  detailed  explanation  of  how 
the  matrix  inversion  subroutine  functions  is  located  in  Appendix  E. 

6.  Error  Identification,  Figure  4.12 

Line  9900  prints  a  message  indicating  that  C2  is  not  invertablc.  Line  9900  is 
based  upon  the  assumption  that  a  division  by  zero  error  in  the  inversion  subroutine 
means  that  C2  is  not  invertable.  Line  9905  prints  the  error  code  and  line  number  of 
other  errors.  See  page  217  of  Reference  1  for  an  explanation  of  error  codes. 
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•INVERT  C2 

FORI 1=1T0NZ: FORI 2=1T0NZ:B1( II, 12)  =C  2(11, 12  )  :NEXTI2 :NEXTI1 

F0RI1=NZ+1T02#NZ: FORI2=lTONZ 

IFI1=I2+NZTHENB1< 12,11 )=1ELSEB1< 12,11 )=0 

NEXTI2:NEXTI1 

FORIl=lTONZ 

ML=1/B1( 11,11 ) : F0RI3=1T02#NZ:B1( II, IS  )=B1( 11,13 )*ML :NEXTI3 
IFI1=NZTHEN9865 

F0RI2=I1+1T0NZ:IFB1( 12,11  )=0THEN9860 
ML=-Bll 12,11) 

F0RI3=I1T02*NZ:B1(  I2',I3  )=B1<  12,13  )  +  (ML#Bl(  11,13  ) ):NEXTI3 

NEXTI2:NEXTI1 

F0RI1=N2T02STEP-1 

F0RI2=I1-1T01STEP-1 :IFB1(  12,11 )=0THEN9885 
ML=-B1( 12,11) 

F0RI3=1T02*NZ:B1( 12,13 )=B1< 12,13  )+<ML*Bl< 11,13 )) :NEXTI3 

NEXTI2 :NEXTI1 

FORI 1 = ITONZ : FORI 2 = 1TONZ 

C2( 12,11  )=B1( 12 ,11+NZ ) :NEXTI2 :NEXTI1 

MI=1: RETURN 


Figure  4.11  Inversion  Subroutine  For  C2. 


9900  IFERL>9800ANDERR=11THENPRINT"! ! TERROR:  C2  Is  Not  Invertabla* • ! " : END 
9905  PRINT"Error  Coda" » ERR >"In  Lina" tERL : END 


Figure  4.12  Error  Identification. 


V.  MODELS  OF  COMBAT  USING  LANCHESTER  EQUATIONS 
A.  GENERAL 

This  program  is  an  example  of  a  time  step  force  attrition  simulation  using 
Lanchester  equations.  The  scenario  is  that  there  are  two  sides,  refered  to  hereafter  as 
the  attacker  and  the  defender.5  The  battle  may  be  broken  into  phases  if  some  model 
parameters  change  during  the  course  of  the  battle.  Each  side  has  a  fixed  number  of 
weapon  types  throughout  the  battle.  The  number  of  attacking  and  defending  weapon 
types  may  be  different.  The  following  characteristics  must  be  specified  for  each 
weapon  type  on  each  side  and  do  not  change  over  the  course  of  the  battle. 

•  The  number  of  weapons  at  the  start  of  the  battle. 

•  The  break  point,  i.e.  the  fraction  of  the  starting  number  of  weapons  which,  if 
reached,  would  cause  the  battle  to  end.  For  example,  if  the  attacker  would 
withdraw  if  50%  of  a  certain  weapon  type  was  lost,  then  the  breakpoint  for  that 
weapon  type  would  hr  j. 

•  The  Lanchester  weapon  characteristic,  i.e.  whether  it  is  a  square  law,  linear  law, 
logarithmic  law  weapon,  or  a  hybrid. 

The  following  characteristics  of  each  weapon  type  may  change  from  phase  to  phase  of 
the  battle. 

•  The  time  required  to  complete  that  phase. 

•  The  rate  at  which  replacements  arrive  for  each  weapon  type. 

•  Attrition  coefficients,  i.e.  the  rate  at  which  a  weapon  type  is  attritted  by  each 
opposing  weapon  type  in  a  particular  battle  phase. 

If  there  are  not  more  than  five  weapon  types  per  side,  the  program  prints  a  dynamic 

display  to  the  M100  screen  showing  for  each  weapon  type  the  fraction  of  starting 

strength  tnat  has  survived  and  the  breakpoint.  Regardless  of  the  number  of  weapon 

types,  the  program  creates  an  output  file  which  shows  the  number  of  survivors  at  the 

end  of  each  phase,  which  weapon  reached  its  breakpoint  first,  and  the  number  of 

survivors  at  the  end  of  the  battle. 


purposes. 


terms  attacker  and  defender  arc  arbitrary'  and  arc  used  only  for  notational 
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B.  REPRESENTING  LANCHESTER  CHARACTERISTICS  OF  EACH  WEAPON 

The  program  includes  provisions  for  calculating  attrition  with  traditional 
Lanchester  square  law  and  linear  law  equations  or  using  a  Helmbold  equation. 
Traditional  Lanchester  attrition  uses  different  functions  for  different  types  of  attrition. 
The  linear  law  functions  (see  Equations  5.1)  for  changes  in  strength  with  respect  to 
time  are  used  to  model  fire  that  is  aimed  at  a  general  area  in  which  targets  are  believed 
to  be  located  An  example  of  a  linear  law  weapon  would  be  artillery  fire  without 
correction  by  an  observer. 


dx/dt  =  -axy,  and  (eqn  5.1) 

dy/dt  =  -byx 

The  square  law  functions  (see  Equations  5.2)  are  used  for  fire  that  is  aimed  at  a  point 
instead  of  a  general  area. 

dx/dt  =  -ay,  and  (eqn  5.2) 

dy/dt  =  -bx 

The  logarithmic  law  functions  (see  Equations  5.3)  are  used  to  model  non-combat  losses 
such  as  disease. 


dx/dt  =  -ax,  and  (eqn  5.3) 

dy/dt  =  -by 

The  reasons  for  using  these  equations  for  modeling  these  types  of  attrition  are 
explained  in  [Ref.  9:Chapter  2].  The  limitation  of  using  Equations  5.1,  5.2,  and  5.3  is 
that  they  restrict  the  model  to  three  discrete  weapons  types.  However,  there  may  be 
weapons  that  do  not  fall  cleanly  into  any  of  the  three  weapon  types.  For  example, 
artillery  fire  that  is  corrected  by  an  observer  may  have  a  mixture  of  linear  law  and 
square  law  characteristics.  The  traditional  Lanchester  equations  also  assume  that  the 
full  fire  power  of  all  the  weapons  on  one  side  can  be  brought  to  bear  against  all  the 
targets  on  the  opposing  side.  However,  in  battles  where  one  side  vastly  outnumbers 
the  other,  or  when  there  are  significant  terrain  masking  or  reaction  time  effects,  this 
assumption  is  not  valid.  To  account  for  these  limitations,  R.  Helmbold  proposed  a  set 
of  modified  Lanchester  equation  in  Reference  8.  Hclmbold's  equations  and  their 
empirical  validation  are  also  discussed  in  [Ref.  9:pp.  174-181  and  Footnote  2.40].  The 
special  cases  of  the  Helmbold  equations  used  in  this  program  are  those  listed  as 
Equations  5.4. 
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(eqn  5.4) 


dx/dt  =  -a(x/y)toy  =  -ax<0x1‘0> 
dy/dt  =  -b(y/x)®)x  =  -by^x1'® 

They  include  an  additional  factor,  (x/y)0*  or  (y/x)10.  If  ct)=  1,  the  result  is  the 
logarithmic  law  equation.  If  (0  =  0,  the  result  is  the  square  law  equation.  If  (0=  .5,  the 
result  is  an  equation  which  behaves  quite  like  the  linear  law  equation,  to  may  be  set  at 
any  value  between  zero  and  one  to  characterize  weapons  which  do  not  fall  neatly  into 
one  of  the  discrete  weapon  types. 

In  the  traditional  Lanchester  linear  law  equation  dx/dt  varies  proportionately  to 
the  number  of  X  survivors  and  the  number  of  Y  survivors.  In  the  Helmbold  equation 
when  (0=  .5,  dx/dt  varies  proportionately  to  the  square  root  of  the  number  of  X  and  Y 
survivors.  If  dx/dt  were  equal  for  the  Lanchester  linear  law  and  the  Helmbold 
equation  with  to  =  .5,  then  a^xy  =  a^(xy)'^  where  a^  and  a j  j  are  the  Lanchester  and 
Helmbold  attrition  coefficients  respectively.  Thus,  a^.j  =  a^(xy)'^.  Therefore,  there  is 
no  single  value  of  a^  that  will  be  equivalent  to  a  value  of  a^  throughout  an  entire 
simulation  because  x  and  y  are  changing.  However,  the  general  shapes  of  the  functions 
kj=  xy  and  k2==(xy)'~’  are  similar  enough  that  they  cause  attrition  to  behave  about 
the  same  in  both  circumstances.  A  plot  of  contour  lines  of  dx/dt  =  1,  2,  and  3  for 
linear  law  and  comparable  Helmbold  equations  is  shown  in  Figure  5.1  where  aj  j  =  .02- 
and  a^=.0002.  The  contour  dx/dt  =  2  is  the  same  for  both  formulations.6  For 
dx/dt  <  2  the  Helmbold  equations  produce  smaller  attrition  rates  than  do  the  traditional 
Lanchester  equations.  For  dx/dt  >2  the  Helmbold  equations  produce  larger  attrition 
rates  than  do  the  traditional  lanchester  equations. 

The  program  includes  BASIC  code  for  both  the  traditional  Lanchester  and 
Helmbold  equations.  To  use  the  traditional  Lanchester  equations,  lines  223-224  and 
233-234  should  be  active,  and  lines  225  and  235  should  be  commented  out  or  deleted. 
To  use  the  Helmbold  equations,  lines  225  and  235  must  be  active,  and  lines  223-224 
and  233-234  must  be  commented  out  or  deleted.  The  weapon  type  parameters  in  the 
input  file  must  match  the  equations  which  are  active  in  the  program. 


6The  Helmbold  coutours  were  jittered  slightly  to  avoid  being  overprinted  at 
dx/dt  =  2  by  the  Lanchester  contour. 


Lanchester— Bold;  Helmbold— Not  Bold 


Figure  5.1  dx/dt  For  Linear  Law  and  Helmbold  Equations. 

EXPLANATION  OF  VARIABLES 

AA(NA,ND)  is  a  matrix  of  the  rates  at  which  attacking  weapons  by  type  are 
attritted  by  each  type  of  defending  weapon  in  a  particular  phase  of  the  battle. 
AB(2,NA)  is  the  matrix  for  the  breakpoints  of  each  attacking  weapon  type. 
Elements  of  the  first  row,  AB(l,i),  are  the  breakpoints  expressed  as  fractions  of 
the  starting  total,  i.e  0  ^  AB(l,i)  ^  1  for  i=  1,2,3...NA.  Elements  of  the  second 
row,  AB(2,i),  are  breakpoints  expressed  as  numbers  of  weapons  of  type  i,  i.e.  0  £ 
AB(2,i)  <  SD(i)  for  i=  1,2,3...NA. 

AT(NA)  is  a  vector  of  tuning  parameters  that  specify  the  Lanchester 
characteristics  of  each  attacking  weapon  type.  If  line  235  is  active  and  lines  233 
and  234  are  commented  out,  then  attrition  is  calculated  using  the  Helmbold 
equation,  see  equation  5.4.  AT(i)  is  to  and  equals  0  for  an  aimed  fire/square  law 
weapon,  .5  for  a  weapon  similar  to  a  Lanchester  area  fire/linear  law  weapon,  and 
1  for  a  source  of  non-combat/logarithmic  law  casualties  to  the  defender.  If  lines 
233  and  234  are  active  and  line  235  is  commented  out,  then  attrition  is  calculated 
using  standard  linear  law  and  square  law  Lanchester  equations,  see  equations  5.1 
and  5.2.  AT(i)  equals  1  for  linear  law  weapons  and  2  for  square  law  weapons. 
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BB(ND,NA)  is  a  matrix  of  the  rates  at  which  one  item  of  each  attacking  weapon 
type  attrits  each  defending  weapon  type  in  a  particular  phase  of  the  battle. 
DB(2,ND)  is  the  matrix  for  the  breakpoints  of  each  defending  weapon  type. 
Elements  of  the  first  row,  DB(l,j),  are  the  breakpoints  expressed  as  fractions  of 
the  starting  total,  i.e  0  ^  DB(l,j)  <  1  for  j=  1,2,3...ND.  Elements  of  the  second 
row,  DB(2,j),  are  breakpoints  expressed  as  numbers  of  weapons  of  type  j,  i.e.  0 
<  DB(2,j)  <  SD(j)  for  j=  1,2,3. ,.ND. 

DT(ND)  is  a  vector  of  tuning  parameters  that  specify  the  Lanchester 
characteristics  of  each  defending  weapon  type.  If  line  225  is  active  a. id  lines  223 
and  224  are  commented  out,  then  attrition  is  calculated  using  the  Helmbold 
equation,  see  equation  5.4.  DT(j)  is  to  and  equals  0  for  an  aimed  fire/square  law 
weapon,  .5  for  a  weapon  similar  to  a  Lanchester  area  fire/linear  law  weapon,  and 
1  for  a  source  of  non-combat/logarithmic  law  casualties  to  the  attacker.  If  lines 
223  and  224  are  active  and  line  225  is  commented  out,  then  attrition  is  calculated 
using  standard  linear  law  and  square  law  Lanchester  equations,  see  Equations  5.1 
and  5.2.  DT(j)  equals  1  for  linear  law  weapons  and  2  for  square  law  weapons. 
11,12,13,  and  14  are  loop  counters. 

MD  is  the  maximum  of  NA  and  ND. 

NA  is  the  number  of  attacking  weapon  types. 

ND  is  the  number  of  defending  weapon  types. 

NI  is  the  number  of  intervals  into  which  a  phase  of  battle  is  broken. 

NP  is  the  number  of  phases  of  battle. 

OA(NA)  is  the  vector  of  horizontal  pixel  positions  for  the  current  screen 
representation  of  the  fraction  of  surviving  attackers  by  weapon  type. 

OD(ND)  is  the  vector  of  horizontal  pixel  positions  for  the  current  screen 
representation  of  the  fraction  of  surviving  defenders  by  weapon  type. 

QA(2,NA)  holds  the  surviving  number  of  each  attacking  weapon  type.  QA(l,i)  is 
the  number  at  the  start  of  a  time  increment,  DT.  QA(2,i)  is  the  number  after 
attrition  by  each  defending  weapon  is  subtracted. 

QD(ND)  holds  the  surviving  number  of  each  defending  weapon  type  after 
attrition  by  each  attacking  weapon  is  subtracted. 

SA(N’A)  is  the  number  of  attacking  weapons  by  type  at  the  start  of  the  battle. 
SD(ND)  is  the  number  of  defending  weapons  by  type  at  the  start  of  the  battle. 


•  TF  is  the  termination  flag.  TF  =  0  means  a  breakpoint  has  not  been  reached. 
TF=  1  means  a  breakpoint  has  been  reached. 

•  TP  is  the  top  pixel  position  for  a  rectangle  or  line  drawn  on  the  graphical  display. 

•  TT  is  the  total  time  that  has  elapsed  in  the  battle  up  to  the  current  time 
increment.  When  a  breakpoint  is  reached,  or  when  all  phases  of  the  battle  are 
completed,  TT  is  time  length  of  the  battle  reported  to  LANOLT.DO. 

D.  INPUT 

All  input  to  the  program  is  entered  into  a  RAM  file,  LANIN. DO.  The  following 
parameters  must  be  entered  in  the  following  order  and  do  not  change  between  battle 
phases. 

•  NP.NA.ND 

•  QA(NA) 

•  QD(ND) 

•  AB(NA) 

•  DB(ND) 

•  AT(NA) 

•  DT(ND) 

The  following  parameters  must  be  entered  in  the  following  order  for  each  phase. 

•  TT,  NT 

•  AR(NA) 

•  DR(NR) 

•  AA(NA.ND) 

•  BB(ND,NA) 

An  example  of  an  input  file  is  shown  in  Figure  5.2 
The  situation  to  be  simulated  using  the  parameters  in  Figure  5.2  is  as  follows.  The 
battle  has  two  phases. 

1.  Parameters  Common  To  Both  Phases 

The  attacker  has  two  weapon  types,  A-,  i  =  1  and  2.  The  defender  has  three 
weapon  types,  D,,  i  =  1,  2,  and  3.  The  attacker  starts  with  200  Ai's  and  100  A->'s. 
The  defender  starts  with  100  Dj's,  200  D2's,  and  100  D^'s.  The  battle  will  end  when 
the  first  attacking  or  defending  weapon  type  is  attritted  to  50%  of  its  initial  strength, 
i.e.  the  breakpoint  for  all  weapon  types  is  .5.  If  this  simulation  is  to  be  run  using 
traditional  Lanchcster  equations,  both  attacking  weapons  are  square  law  weapons.  Dj 
is  an  linear  law  weapon,  and  D2  and  D-j  are  square  law  weapons. 


Figure  5.2  Sample  Input  File. 


2.  Situation  For  Phase  One 


The  length  of  the  first  phase  is  five  hours.  The  program  will  break  that  period 
into  50  segments  for  computational  purposes.  Replacements  arrive  for  both  attacking 
weapon  types  at  an  average  rate  of  one  weapon  per  hour.  The  replacement  rate  for  all 
three  defending  weapon  types  averages  .8  weapons  per  hour.  Aj's  are  attritted  by  each 
Dj  at  a  rate  of  .00012  per  hour  per  surviving  Al.  Aj's  are  attritted  by  each  D2  at  a 
rate  of  .014  per  hour.  Rates  at  which  the  remaining  Aj  are  attritted  by  each  Dj  are 
listed  through  the  end  of  the  next  line.  Dj's  are  attritted  by  each  Aj  at  a  rate  of  .019 
per  hour.  Dj's  are  attritted  by  each  A2  at  a  rate  of  .017  per  hour.  Rates  at  which  the 
remaining  Dj  are  attritted  by  each  Aj  are  listed  through  the  end  of  the  next  two  lines. 

3.  Situation  For  Phase  Two 

The  length  of  the  second  phase  is  ten  hours.  The  program  will  break  that 
period  into  50  segments  for  computational  purposes.  Replacements  arrive  for  both 
attacking  weapon  types  at  the  rate  of  .5  weapons  per  hour.  The  replacement  rate  for 
all  three  defending  weapon  types  averages  .4  weapons  per  hour.  Aj's  are  attritted  by 
each  Dj  at  a  rate  of  .00018  per  hour  per  surviving  Al.  Aj's  are  attritted  by  each  D2 
at  a  rate  of  .021  per  hour.  Rates  at  which  the  remaining  Aj  are  attritted  by  each  Dj 


are  listed  through  the  end  of  the  next  line.  Dj's  are  attritted  by  each  Aj  at  a  rate  of 
.0285  per  hour.  Dj's  are  attritted  by  each  A2  at  a  rate  of  .0255  per  hour.  Rates  at 
which  the  remaining  D;  are  attritted  bv  each  A;  are  listed  through  the  end  of  the  next 
two  lines. 

E.  OUTPUT 

The  program  produces  two  types  of  output: 

•  An  output  file.  LANOUT.DO,  which  lists  the  status  of  each  weapon  type  at  the 
end  of  each  phase  and  at  the  end  of  the  battle  and  which  weapon  types  have 
gone  below  their  breakpoints. 

•  A  dvnamic  graphical  display  to  the  screen  of  the  M100  which  shows  the  fraction 
of  the  starting  strength  of  each  weapon  type  which  has  survived  until  that  time 
interval  in  the  simulation. 

The  graphical  display  consists  of  a  rectangle  on  the  M100  screen  for  each 
attacking  and  defending  weapon  type.  Each  rectangle  is  100  pixels  wide  and  five  pixels 
high.  Each  pixel  in  the  horizontal  direction  represents  one  percent  of  the  starting 
strength  of  the  weapon  type  represented  by  that  particular  rectangle.  In  each  rectangle 
is  a  vertical  line  showing  the  breakpoint  for  that  weapon  type  as  a  fraction  of  starting 
strength.  As  the  simulation  progresses,  another  vertical  line  in  each  rectangle  is 
updated  showing  the  fraction  of  survivors  for  that  weapon  type.  The  rectangles  are 
arranged  in  two  columns,  one  for  attacking  and  one  for  defending  weapon  types. 
Weapon  type  numbers  are  printed  to  the  left  of  the  rectangles.  If  the  replacement  rate 
drives  the  number  of  survivors  over  the  starting  strength  for  some  weapon  type,  the 
vertical  line  indicating  the  fraction  of  survivors  will  stay  at  the  100%  level  and  an 
asterisk  will  be  printed  to  the  left  of  the  corresponding  rectangle.  In  addition  to  the 
rectangles  there  is  a  printed  line  at  the  bottom  of  the  screen  which  tells  the  operator  on 
what  phase  and  time  interval  the  simulation  is  currently  working. 

F.  EXPLANATION  OF  PROGRAM  COMPONENTS 
A  complete  listing  of  this  program  is  at  Appendix  C. 

1.  Initialization  Section,  Figure  5.3 

Line  120  sets  the  number  of  files  to  two  and  opens  the  input  file,  LANIN. DO. 
Line  121  opens  the  output  file,  LANOUT.DO  and  enters  the  number  of  phases  in  the 
battle,  NP,  and  the  number  of  attacking  and  defending  weapon  types,  NA  and  ND. 
Line  122  defines  MD  as  the  maximum  of  NA  and  ND.  If  both  sides  have  five  or  fewer 
weapons  types,  line  124  sets  SF  =  1,  indicating  that  the  screen  of  the  Ml 00  is  big 
enough  to  handle  the  graphical  display  generated  during  the  simulation.  If  cither  side 
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Figure  5.3  Initialization  Section. 

has  more  than  five  weapon  types,  the  graphical  display  will  not  be  generated.  Lines 
130-131  dimension  the  matrices  used  in  the  program. 

2.  Entering  Common  Parameters,  Figure  5.4 


132  'Enter  Initial  Quantities  of  Wpns>  Break  Points  And  Wpn  Types. 

134  FORI 2=1T0NA : INPUT#1 ,QA<  2,12 ) :SA( 12 )=QA(  2,12 ) :OA(  12 )=127:NEXTI2 

135  FORI 2=1TONO : INPUT# 1 ,QD < 1 2  ) : SO < 1 2  )=QD ( 1 2  ) : OD 1 12  ) =238 : NEXTI 2 

136  F0RI2=1T0NA:INPUT#1,AB< 1,12 ): ABC  2,12  )=AB( 1,12  )*QA( 2,12 ) :NEXTI2 

137  F0RI2=1T0ND :INPUT#1,D8( 1,12  ):DB( 2,12  )=DB( 1,12  )*QD( 12 ) : NEXTI 2 

138  F0RI2=1T0NA:INPUT#1,AT( 12  )  :NEXTI2 : F0RI2=1T0ND:INPUT#1,DTI 12  ): NEXTI 2 
140  TM=0 : IFSF =1THENGQSUB600 


Figure  5.4  Entering  Common  Parameters. 

This  section  enters  the  rest  of  the  parameters  that  do  not  change  from  phase 
to  phase  of  the  battle  and  initializes  the  variables  controlling  the  graphical  display. 
Line  134  puts  the  starting  quantity  of  Aj  into  QA(2,i)  and  SA(i)  and  sets  OA(i)  to  127 
which  indicates  that  100%  of  Aj  are  surviving  at  the  start  of  the  simulation.  Line  135 
performs  the  same  function  for  Dj  that  line  134  performed  for  Aj.  Lines  136  and  137 
enter  the  fractional  breakpoints  for  Aj  and  Dj  respectively  into  AB(l,i)  and  DB(  1  ,j). 
Lines  136  and  137  also  compute  the  breakpoints  in  terms  of  numbers  of  surviving 
weapons,  placing  them  in  AB(2,i)  and  DB(2,j).  Line  138  enters  the  Lanchcstcr 
characteristic  parameters,  AT(i)  and  DT(j),  for  each  weapon  type.  Line  140  sets  TM, 
which  keeps  track  of  the  time  until  the  end  of  the  battle,  to  zero.  If  the  graphical 
display  to  the  screen  is  to  be  used,  line  140  also  calls  subroutine  600  which  sets  up  the 
output  screen. 


3.  Initialization  For  Each  Phase,  Figure  5.5 
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143  FORIl=lTONP : PRINTR2 » "STARTING  PHASE**»I1 

145  'Enter  Timo  Spent  In  Phase  II  and  #  of  Intervals 

146  INPUT81 ,TT  >NI :DT=TT/NI 

150  'Enter  Replacement  Rates  And  Attrition  Coefficient  Matrices 
152  F0RI2=1T0NA: INPUTR1 , AR1 12  )  :NEXTI2 : F0RI2  =  1T0ND : INPUTRl ,0R( 12  )  :NEXTI2 
154  F0RI2=1T0NA : FORI3=1TONO : INPUTR1 , AA( 12,13  )  :NEXTI3 :NEXTI2 
158  F0RI2=1T0ND : F0RI3=1T0NA : INPUTS1 ,BB( 12,13  ) :NEXTI3 :NEXTI2 


Figure  5.5  Initialization  For  Each  Phase. 

Line  143  sets  II,  the  phase  counter,  and  prints  a  heading  to  the  output  file. 
Line  146  enters  the  length  of  the  current  phase,  TT,  and  the  number  of  intervals,  NI, 
into  which  the  phase  will  be  broken  and  computes  the  length  of  each  interval,  DT. 
The  choice  of  NT  is  a  compromise  between  two  competing  objectives:  accuracy  and 
time  required  to  complete  the  simulation.  As  NT  becomes  larger,  DT  becomes  smaller 
and  the  simulation  more  closely  approximates  the  continuous,  mutual  attrition  that  is 
the  basis  of  Lanchester  equation  theory.  If  NI  is  small  and  DT  is  large,  then  the 
simulation  tends  to  discount  the  attrition  which  takes  place  during  a  phase  because  the 
program  assumes  force  levels  are  constant  throughout  an  interval,  DT.  However,  if  NI 
is  to  large,  the  time  required  to  run  the  simulation  increases  linearly.  The  relationship 
between  accuracy  and  speed  is  also  a  function  of  the  attrition  coefficients  and  number 
of  weapon  types  on  each  side. 

Line  152  enters  the  replacement  rates  for  each  weapon  type.  Lines  154-156 
enter  the  attrition  coefficients  for  each  weapon  type. 

4.  Attrition  Calculations  For  A  Phase,  Figure  5.6 

This  section  breaks  a  phase  into  NT  intervals  of  length  DT,  calculates  the 
attrition  during  each  interval,  and  tests  whether  that  attrition  has  caused  some  weapon 
type  to  reach  its  breakpoint. 

Line  202  sets  the  interval  counter,  12,  adds  the  length  of  the  interval  to  the 
length  of  the  battle,  TM,  and  prints  a  message  to  the  screen  telling  the  operator  what 
phase  and  interval  is  currently  being  processed. 

Lines  222-227  calculate  the  attrition  to  attackers  based  upon  the  quantities  of 
each  weapon  surviving  at  the  start  of  the  interval.  QA(2,i)  holds  the  current  quantity 
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200  'Fight  Phase  II. 

202  F0RI2=1T0NI :TM=TM+0T : PRINT3241 » "Phase >11 i" ,  Increment" tI2>"out  of'lNI 
210  'Fight  Time  Increment  0T. 

220  ‘Update  number  of  attackers 

222  F0RI3=1T0NA:QAI1,I3)=QAI2,I3):NEXTI3:F0RI3=1T0NA:F0RI4=1T0ND 

223  I*I>T(I4)=1THENQA<  2,13  )=QA<  2,13  )-AA(  13,14 )*QD<  14 >*QA(  2,13  )*DT:G0T0226 

224  QA<  2,13 )=QA<  2,13 )-AA( 13 ,14  )*QD< 14  )*DT 

225  *QA(  2,13  )=QA(  2,13  )-AAl  13 ,14  )*( QAt  2,13  )/QD<  14  )  )ADT<  14  )«QD<  14  )*dt 

226  NEXTI4:QA(  2,13  )=«A<  2,13  J+ARI  13  )*0T: IFSF=1THENGOSUB650 

227  NEXTI3 

230  'Update  number  of  defenders 

232  F0RI3=1T0ND:F0RI4=1T0NA 

233  IFATI 14  )=1THENQD( 13  )=QD< 13  )-BB( 13,14  >*QD( 13  )*QAt 1,14 >*DT:G0T0236 

234  QDI  13 )=QD< 13 )-BBt 13 ,14 )#QA( 1,14  )*DT 

235  'QD(  13  )=QD(  13  )-B8(  13,14  )*(  Q01 13  )/QA(  1,14  )  )A  ATI  14  )*QA(  1,14  )*dt 

236  NEXTI4:QD( 13  )=QD( 13  )+0R( 13  )*OT:IFSF=1THENGOSUB660 

237  NEXTI3 

240  GOSUB300:NEXTI2 

242  IFIl=NPTHENGOSUB350:CLS:PRINT"Output  is  in  file  LAN0UT.D0.":END 
245  PRINTP2, "Status  After  Phase"lIl:GOSUB361:NEXTIl 


Figure  5.6  Attrition  Calculations  for  a  Phase. 

of  Aj  and  is  therefore  decremented  by  lines  222-227.  Attrition  to  each  Dj  should,  for 
consistency  with  the  attrition  to  the  A-'s,  also  be  based  upon  the  quantity  of  Aj's 
surviving  at  the  beginning  of  the  interval.  Therefore,  the  quantity  of  each  A-  surviving 
at  the  beginning  of  the  interval  is  saved  by  line  222  in  QA(l,i)  for  use  during  the 
attrition  calculations  for  Dj.  Line  222  also  sets  the  Aj  counter,  13,  and  the  Dj  counter, 
14. 

If  the  simulation  is  to  be  done  with  traditional  Lanchester  linear  law  and 
square  law  equations,  lines  223-224  must  be  active  and  line  225  must  be  commented 
out  or  deleted.  If  the  simulation  is  to  be  done  with  a  Helmbold  equation,  lines  223-224 
must  be  commented  out  or  deleted  and  line  225  must  be  active.  The  program 
displayed  in  Figure  5.6  has  the  Helmbold  equations  commented  out.  Line  222 
calculates  the  attrition  of  A^  by  Dj  based  upon  a  linear  law,  Equation  5.1.  Line  223 
calculates  attrition  based  upon  the  square  law,  Equation  5.2.  Whether  the  linear  law 
or  square  law  is  used  is  based  upon  AT(i)  and  is  determined  in  line  222.  If  line  225 
were  active,  it  would  calculate  attrition  using  the  Helmbold  equation,  Equation  5.4. 
After  attrition  by  each  Dj  is  calculated,  line  226  adds  the  replacements  for  Aj  received 
during  the  interval.  Line  226  also  calls  subroutine  650  which  updates  the  graphical 
display  to  reflect  the  attrition  to  each  Aj. 


prjFH 


Lines  232-237  calculate  the  attrition  to  defenders  based  upon  the  quantities  of 
each  weapon  surviving  at  the  start  of  the  interval.  QD(i)  holds  the  current  quantity  of 
Dj.  Line  232  sets  the  Dj  counter,  13,  and  the  Aj  counter,  14.  If  the  simulation  is  to  be 
done  with  traditional  Lanchester  linear  law  and  square  law  equations,  lines  233-234 
must  be  active  and  line  235  must  be  commented  out  or  deleted.  If  the  simulation  is  to 
be  done  with  a  Helmbold  equation,  lines  233-234  must  be  commented  out  or  deleted 
and  line  235  must  be  active.  The  program  displayed  in  Figure  5.6  has  the  Helmbold 
equations  commented  out.  Line  232  calculates  the  attrition  of  Dj  by  Aj  based  upon  a 
linear  law,  Equation  5.1.  Line  233  calculates  attrition  based  upon  the  square  law, 
Equation  5.2.  Whether  the  linear  law  or  square  law  is  used  is  based  upon  DT(i)  and  is 
determined  in  line  232.  If  line  235  were  active,  it  would  calculate  attrition  using  the 
Helmbold  equation,  Equation  5.4.  After  attrition  by  each  Aj  is  calculated,  line  236 
adds  the  replacements  for  Dj  received  during  the  interval.  Line  236  also  calls 
subroutine  660  which  updates  the  graphical  display  to  reflect  the  attrition  to  each  Dj. 
Line  240  calls  subroutine  300  to  check  for  whether  a  breakpoint  was  reached  during 
the  interval.  If  no  breakpoint  was  reached,  a  new  interval  is  begun. 

If  all  the  intervals  in  the  last  phase  are  completed  without  reaching  a 
breakpoint,  line  242  calls  subroutine  350  which  prints  the  status  at  the  end  of  the  battle 
to  the  output  file.  If  the  current  phase  is  not  the  last  phase,  then  line  245  prints  a 
header  to  the  output  file,  calls  subroutine  361  which  prints  the  status  at  the  end  of  the 
current  phase  to  the  output  file,  and  starts  the  next  phase. 

5.  Breakpoint  Subroutine,  Figure  5.7 

This  section  determines  whether  any  weapons  have  reached  their  breakpoints, 
and,  if  so,  prints  that  information  in  the  output  file,  and  ends  the  program.  Line  320 
sets  TF  =  0,  indicating  that  no  breakpoints  have  been  reached.  Line  320  then  starts  a 
loop  which  tests  whether  any  Aj  have  reached  their  breakpoints.  If  so,  then  lines 
322-324  print  a  message  to  the  output  file  specifying  the  weapon  type,  the  breakpoint, 
and  the  quantity  of  that  weapon  type  that  survived.  Lines  335-340  test  whether  any  Dj 
have  reached  their  breakpoints,  and  if  so,  print  a  message  to  that  effect  to  the  output 
file.  If  no  breakpoints  have  been  reached,  then  line  340  returns  control  to  the  main 
program  to  begin  attrition  calculations  in  the  next  time  interval. 

The  default  battle  termination  criterion  is  that  at  least  one  weapon  type  must 
be  below  its  individual  breakpoint  at  the  end  of  a  time  interval.  However,  the  operator 
may  wish  to  edit  the  program  before  running  it,  adding  more  sophisticated  termination 


300  'Check  Whether  Breakpoint  is  reached. 

320  TF=0:FORI3=1TONA:IFQA< 2,13  )>AB<  2,13 JTHEN325 

322  TF=1:PRINT#2, "Attacker  Hpn">I3>“Is  Below  Breakpoint" 

323  PRINT#2,"  Bp  =“>  :PRINT#2, USING"####. ##">AB(  2,13  )> 

324  PRINT tt2 ,"  Current  Level  =“ >:PRINT#2, USING"####. ##">QA< 2,13 ) 

325  NEXT I 3 

335  FORI3=1TOND : IFQD( 13 )>DB<  2,13 ITHEN340 

337  TF =1 : PRINTP2 , "Defender  Hpn">I3»"Is  Below  Breakpoint" 

338  PRINT tf2,"  Bp  ="  J  :PRINT#2, USING"####.  ##"  >DB<  2,15  )> 

339  PRINT#2 »"  Currept  Level  =" j :PRINT#2, USING"####. ##" >QD( 13 ) 

340  NEXTI3:IFTF=0THENRETURN 

350  PRINT#2, PRINT#2,"":PRINT#2, "SUMMARY  AT  END  OF  BATTLE" 

351  PRINT# 2 , : PRINT#2 ,"T ime  Elapsed  During  Battle  ="> 

352  PRINT#2 , USING”#### . ##" } TM : PRINT#2 , : GOSUB361 
355  CLS : PRINT"0utput  is  in  -file  LANOUT  .  DO" : END 

361  PRINT#2,"  Att  Wpn  Breakpoint  Current  Level" 

363  F0RI3=1T0NA : PRINT#2 , USING"######" >13 v  . 

364  PRINT#2, USING"##########. ##">AB( 2,13 )>QA< 2, 13 > :NEXTI3 :PRINT#2,"" 

366  PRINT#2,“  Oef  Wpn  Breakpoint  Current  Level" 

367  F0RI3=1T0ND :  PRINT#2  .USING'1######"  >13  > 

368  PRINT#2, USING"##########. ##" >DBI 2 ,13 ) >GD(  13 ):NEXTI3:PRINT#2,"" : RETURN 


Figure  5.7  Subroutine  To  Test  For  Breakpoints. 


criteria  to  the  default  criteria.  For  example,  if  the  operator  wants  the  battle  to 
terminate  when  Aj  or  A2  reach  half  their  starting  strength  or  when  Aj  reaches  60% 
and  A2  reaches  70%  of  their  starting  strength  the  operator  should: 

•  Put  .5  as  the  individual  breakpoints  for  Aj  and  A2  in  the  input  file  and 

•  Put  the  program  lines  shown  in  Figure  5.8  into  the  program  after  line  325. 


330  IFQAI  2,1 )/SA(  1 )>. 60RQAI 2,2 )/SAl 2  )>. 7THEN335 

331  TF=1:PRINT#2, "Special  Termination  Criterion  Met" 


Figure  5.8  Example  Of  An  Additional  Termination  Criterion. 


If  a  breakpoint  has  been  reached,  then  lines  350-368  print  the  status  of  both 
sides  at  the  end  of  the  battle.  That  end  of  battle  status  report  includes  a  header,  time 
elapsed  during  the  battle,  and  a  list  of  attacking  and  defending  weapon  types  with  their 
breakpoints  and  number  of  survivors.  Lines  361-368  are  called  as  a  subroutine  from 
line  352  because  lines  361-368  are  also  used  to  print  the  summary  at  the  end  of  each 
phase  and  can  therefore  not  terminate  the  program. 


6.  Graphics  Display  Initialization  Subroutine,  Figure  5.9 


600  'Set  up  output  screen 

610  PRINT"Hpn  #  Attacker  Defender" 

620  F0RI1=1T0MD: PRINTUSING"##" >11 

623  TP=2+I1*8 

625  IFI1>NATHEN630 

627  LINE! 18, TP  )-( 119.TP+4  > ,1,B:BP=18+INT( 100«AB( 1,11 ) ) 

628  LINE1 BP-1 ,TP+1 )-( BP ,TP»3  ) , 1 ,B 
630  IFI1>N0THEN635 

632  LINE) 138, TP )-( 239.TP+4 ) ,1 ,B:BP=138+INT( 100«DB( 1,11 ) ) 

633  LINE! BP-1, TP+l)-< BP, TP+3),1,B 
635  NEXTI1: RETURN 


Figure  5.9  Graphics  Display  Initialization  Subroutine. 

The  graphics  display  consists  of  a  rectangle  on  the  M100  screen  for  each 
attacking  and  defending  weapon  type.  Each  rectangle  is  100  pixels  wide  and  five  pixels 
high.  Each  pixel  in  the  horizontal  direction  represents  one  percent  of  the  starting 
strength  of  the  weapon  type  represented  by  a  particular  rectangle.  The  rectangles  are 
arranged  in  two  columns,  one  for  attacking  and  one  for  defending  weapon  types.  This 
subroutine  draws  the  rectangles,  labels  the  columns  "Attacker"  or  "Defender",  puts  a 
vertical  line  in  each  box  at  the  breakpoint  for  that  weapon  type,  and  labels  the  rows  of 
boxes  with  the  weapon  type  number. 

Line  610  prints  the  header  on  line  one  of  the  screen.  Line  620  starts  a  loop 
that  writes  the  weapon  type  number,  II,  and  prints  the  rectangles;  Line  623  calculates 
the  vertical  pixel  position,  TP,  for  the  top  of  the  rectangles  of  weapon  type  II.  Line 
625  tes' s  w'hether  to  draw  a  rectangle  next  to  weapon  type  number  1 1  in  the  "Attacker" 
column.  If  so,  the  first  statement  on  line  627  draws  the  rectangle.  The  second 
statement  on  line  627  calculates  the  horizontal  pixel  location  in  the  rectangle  of  the 
breakpoint  for  that  weapon  type.  Line  628  draws  a  double  line  at  the  breakpoint  in 
the  rectangle.  Lines  630-635  test  whether  a  rectangle  should  be  drawn  in  the  defender 
column.  If  so,  they  draw  the  rectangle  and  insert  the  breakpoint  in  the  same  manner 
as  was  done  in  lines  620-628.  Line  635  returns  control  to  the  main  program  when 
there  are  no  more  rectangles  to  be  drawn. 


7.  Updating  The  Graphical  Display,  Figure  5.10 


650  'Update  screen  output  of  attackers 
653  TP=3+I3*8 

655  LINE) OA( 13  )  ,TP  )-( OA< 13  )  ,TP+2  )  ,0 

656  0A( 13  )=18+INT< 100*QAl 2,13  )/SA( 13  ) ) 

657  XFOAI 13 )>118THEN0A< 13 )=X18: PRINTS! I3*40+2 G0T0659 

658  PRINT3I3*40*2,"  " 

659  LINE ( 0A( 13  )  ,TP )-( OA( 13 ) ,TP+3 ) , 1 : RETURN 

660  ‘update  screen  output  of  defenders 
663  TP=3+I3*8 

665  LINE! 0D< 13  )  ,TP  )-( 0D( 13  )  ,TP+2  )  ,0 

666  00< 13  )=138+INT( 100*QD< 13 )/SD(  13 ) ) 

667  IFO0<I3)>238THENOAlI3)=238:PRINT»I3*40+22,,'»,,:GOTO669 

668  PRINT3I3«40+22,“  “ 

669  LINE ( 00 ( 13  )  ,TP )-< 0D( 13 ) ,TP+3 ) ,1: RETURN. 


Figure  5.10  Subroutines  To  Update  The  Graphical  Display. 

This  section  includes  two  subroutines  which  update  the  vertical  line  in  each 
rectangle  which  indicates  the  fraction  of  survivors  for  that  weapon  type.  The 
subroutine  in  lines  650-659  updates  the  attacker  rectangles;  lines  660-669  perform  the 
same  function  for  defender  rectangles.  Line  653  sets  TP,  the  location  of  the  top  pixel 
of  the  vertical  line  for  Aj-j.  Line  655  erases  the  old  vertical  line,  the  horizontal  pixel 
position  for  which  was  stored  in  OA(I3).  Line  656  calculates  the  horizontal  pixel 
position  for  the  new  vertical  line  based  upon  the  fraction  of  the  starting  strength  of 
Aj3's  which  currently  survives.  If  the  reinforcement  rate  exceeds  the  attrition  rate  and 
drives  the  number  of  survivors  over  the  starting  strength  for  Ajj  line  657  holds 
horizontal  position  of  the  vertical  line  at  the  100%  level  and  prints  an  asterisk  next  to 
the  corresponding  rectangle.  If  the  number  of  survivors  is  less  than  the  starting 
strength,  line  658  prints  a  blank  space  next  to  the  corresponding  rectangle.  Line  659 
writes  the  new  vertical  line  to  the  screen  showing  the  fraction  of  Aj^'s  which  survive. 
Lines  660-669  perform  the  same  function  for  Dj  that  lines  650-659  perform  for  Aj. 

G.  EXAMPLE  SIMULATIONS 
1 .  Example  #  1 

The  first  example  uses  the  input  file  shown  in  Figure  5.2  The  output  file  for 
that  simulation  is  shown  in  Figure  5.11. 
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STARTING 

PHASE  1 

Status  Aftar  Phase  1 

Att 

Breakpoint 

Current  Level 

X 

100.00 

173.93 

2 

50.00 

68.55 

Dat  Hpn 

Breakpoint 

Current  Level 

1 

50.00 

79.12 

2 

100.00 

184.53 

3 

50.00 

89.95 

STARTING 

PHASE  2 

Attacker  Hpn  2  Is  Below  Breakpoint 

Bp  = 

50.00  Current 

Level  =  48.85 

SUMMARY  AT  END  OF  BATTLE 

Time  Elapsed  During  Battle  =  7.20  1 

Att  Hpn 

Breakpoint 

Current  Level 

1 

100.00 

157.29 

2 

50.00 

48.85 

Oaf  Wpn 

Breakpoint 

Current  Level 

1 

50.00 

66.24 

2 

100.00 

174.64 

3 

50.00 

84.25 

Figure  5.11  Output  File,  LANOUT.DO,  For  Example  #1. 


2.  Comparing  The  Lanchester  and  Helmbold  Linear  Law  Equations 

Examples  2  and  3  compare  the  differences  between  using  a  traditional 
Lanchester  linear  law  equation  (Example  2)  and  a  Helmbold  equation  (Example  3). 
The  scenereos  for  Examples  2  and  3  share  the  following  elements. 


The  battle  has  only  one  phase  with  a  maximum  length  of  10  which  is  broken  into 
100  intervals.  & 


•  Both  sides  have  two  weapon  types.  Each  weapon  type  has  a  starting  strength  of 


•  The  breakpoints  for  all  weapon  types  are  50%. 

•  All  weapon  types  are  linear  law  weapons.  Attrition  is  calculated  using  Equation 
5.1  for  Example  2  and  using  Equation  5.4  (co=  .5)  for  Example  3. 

•  There  are  no  replacements  for  any  weapon  type. 

The  only  differences  in  the  scenereos  for  Examples  2  and  3  are  the  attrition 
coefficients. 

a.  Example  #2 

The  input  and  output  files,  LANIN. DO  and  LANOUT.DO,  for  Example  2 
are  in  Figure  5.12.  Since  these  attrition  rates  are  for  the  Lanchester  linear  law,  the 
dimensionality  of  the  rates  for  the  attacker  are  (number  of  attacker  casualties)  per 


(number  of  attackers)  per  (number  of  defenders)  per  (unit  time).  The  dimensionality  of 
the  rates  for  the  defender  are  the  same  with  the  rolls  reversed. 


Input  File: 

Output  File: 

12  2 

STARTING  PHASE  1 

100  100 

Attacker  Wpn  1  Is  Below  Breakpoint 

10 0  100 

Bp  =  50.00  Current 

Level  =  49.40 

.5  .5 

.5  .5 

SUMMARY  AT  END  OF  BATTLE 

1  1 

Time  Elapsed  During  Battle  =  5.20 

1  1 

10  100 

Att  Wpn  Breakpoint 

Current  Level 

0  0 

1  50.00 

49.40 

0  0 

2  .50.00 

62.53 

.00075  .00075 

.0005  .0005 

Def  Wpn  Breakpoint 

Current  Level 

.00025  .00025 

1  50.00 

82.17 

.00025  .00025 

2  50.00 

82.17 

Figure  5.12  Input  And  Output  Files  For  Example  #2. 


b.  Example  # 3 

The  input  and  output  files,  LANIN. DO  and  LANOUT.DO,  for  Example  3 
are  in  Figure  5.13.  The  attrition  rates  in  this  example  are  for  the  Helmbold  Equation 
with  to  =  .5,  the  Helmbold  equivalent  of  the  Lanchester  linear  dimensionality  of  the 
rates  for  the  attacker  are  (number  of  attacker  casualties)  per  (attacker)"'’  per 
(defender)'^  per  (unit  time).  The  dimensionality  of  the  rates  for  the  defender  are  the 
same  with  the  rolls  reversed. 

To  generate  Helmbold  coefficients  that  are  comparable  to  the  Lanchester 
linear  law  coefficients,  the  Lanchester  coefficients  must  be  adjusted  by  the  difference  in 
the  dimensionality,  i.e.  multiplied  by  [(number  of  attackers)(number  of  defenders)] 

In  this  example  it  means  multiplying  the  Lanchester  coefficients  by  100. 

The  results  of  Examples  2  and  3  show  good  agreement. 

•  The  simulation  using  the  Helmbold  equations  ended  about  21%  faster  than  did 
the  battle  using  Lanchester  equations.  The  same  weapon  type  reached  its 
breakpoint  first  in  both  cases. 

•  The  differences  between  the  number  of  survivors  for  other  weapon  types  was 
quite  small. 


Input  File: 


Output  File: 


12  2 
100  100 
100  100 
.5  .5 
.5  .5 
.5  .5 
.5  .5 
10  100 
0  0 
0  0 

.0000075 
.000005 
.0000025 
. 0000025 


STARTING  PHASE  1 

Attacker  Wpn  1  Is  Below  Breakpoint 
Bp  =  50.00  Current  Level  =  49.85 

SUMMARY  AT  END  OF  BATTLE 

Time  Elapsed  During  Battle  =  4.10 


Att  Wpn 

Breakpoint 

Current  Level 

1 

50.00 

49.85 

.0000075 

2 

50.00 

64.67 

000005 

Def  Wpn 

Breakpoint 

Current  Level 

.0000025 

1 

50.00 

82.79 

.0000025 

2 

50.00 

82.79 

Figure  5.13  Input  And  Output  Files  For  Example  #3. 


VI.  GEOMETRIC  PROGRAMMING 


A.  GENERAL 

The  program  described  in  this  chapter  solves  nonlinear  programming  problems  in 
which: 

•  The  objective  function  is  to  be  minimized. 

•  The  objective  function  and  constraints  are  posynomials. 

•  The  number  of  terms,7  T,  minus  the  number  of  variables,  N,  must  equal  one. 

•  The  coefficients,8  cm  t,  of  all  terms  must  be  strictly  positive. 

•  All  components  of  the  vector  of  decision  variables,  x>  must  be  strictly  positive  at 
optimality. 

•  Constraints  must  have  the  form  of  a  posynomial  on  the  left  hand  side  that  is  less 
than  or  equal  to  one. 

Geometric  programming  has  the  distinctive  feature  of  calculating  the  optimal 
value  of  the  objective  function  before  the  optimal  values  of  the  decision  variables  are 

calculated.  Geometric  programming  also  produces  weights,  5t,  t=  1,2,3 . T,  associated 

with  each  term.  For  example,  in  applications  where  the  ct  are  prices  and  the  objective 


7The  number  of  terms  in  the  objective  function  and  in  the  constraints. 

8Two  subscripting  systems  are  used  throughout  this  chapter.  The  first  uses  the 
letters  m  and  t  where  t  =  l,2,3,...,Tm  is  the  number  of  the  term  in  the  mth  posynomial. 
m  =  0  refers  to  the  objective  function;  m=  1,2,...  refers  to  the  constraint  numbers.  Tm 
is  the  number  of  terms  in  the  mth  constK.nt,  When  problems  of  only  one  posynomial 
are  being  discussed,  the  m  is  omitted.  The  first  system  also  includes  the  letter  n  to 
identify  components  of  the  decision  variable  vector,  x.  where  n=  1,2,. ...N.  The  second 
system  numbers  the  terms  without  starting  again  at  1  at  the  beginning  of  each 
constraint.  Each  term  is  numbered  t',  t'=  1,2,...,T'  where  T  is  the  number  of  terms  in 
the  objective  function  and  the  constraints.  The  first  term  of  the  objective  function  is 
denoted  by  f  =  1.  The  other  terms  in  the  objective  function  are  then  numbered  from 
left  to  right.  Then  the  terms  in  each  constraint  in  turn  arc  numbered  from  left  to  right. 
For  example,  in  Figure  6.1,  t=  1,  m  =  0  (or  f  =  1)  refers  to  4<)xj\->  and  t=2,  m=  1  (or 
t'  =  5)  refers  to  ^x-flx-j'*'  7. 


function  minimizes  total  cost,  the  5^  for  objective  function  terms  are  the  proportion  of 
cost  that  term  t  contributes  to  optimal  total  cost,  f(x  )•  These  weights  are  invariant 
with  respect  to  the  prices,  ct,  associated  with  each  term. 

1.  Definition  Of  A  Posynomial  Function 

The  function  f|x)  is  posynomial  if  it  has  the  form 

flX)  =  Z  ct  Pt(X)  where  Pt(X)  =  n  xn  n,t  (eqn  6’^ 

t=l  n=l 

where 

•  T  is  the  number  of  terms. 

•  ct  are  positive  scalar  constants. 

•  x  is  the  vector  of  decision  variables,  (xj^.—.x^). 

•  The  only  restriction  on  the  exponents,  am  n^,  is  that  they  be  real  numbers. 

A  posynomial  differs  from  a  polynomial  in  that  the  coefficients  of  a  posynomial  must 
be  strictly  positive  and  its  exponents,  an  t,  need  not  be  positive  integers. 

An  example  of  a  problem  meeting  these  conditions  is  in  Figure  6.1. 


Figure  6.1  Geometric  Programming  Problem  In  Standard  Form. 

Geometric  programming  solves  a  problem  of  this  kind  by  solving  its  dual. 
When,  as  specified  above,  the  number  of  terms  minus  the  number  of  variables  equals 
one,  then  the  problem  has  a  unique  solution.  T  -  N  -  1  is  by  convention  called  the 
degree  of  difficulty.  If  the  degree  of  difficulty  is  greater  than  zero,  then  another 
nonlinear  program  must  be  solved  to  find  the  optimal  5^  .  While  this  new  nonlinear 
program  is  frequently  easier  to  solve  than  the  original  problem,  its  solution  is  beyond 
the  scope  of  this  chapter  which  is  limited  to  problems  with  a  degree  of  difficulty  of 


zero. 


B.  MATHEMATICAL  BASIS  FOR  GEOMETRIC  PROGRAMMING 

The  mathematical  basis  for  geometric  programming  is  summarized  in  [Ref.  10:pp. 
494-522].  and  explained  in  detail  in  Reference  11.  The  following  explanation  is 
provided  for  tutorial  purposes  and  is  an  adaptation  of  the  explanation  in  [Ref.  10:pp. 
496-502].  Notation  in  this  chapter  is  consistent  with  that  used  in  Reference  10. 
a.  Unconstrained  Minimization  Of  Posynomial  Functions 

Historically,  geometric  programming  has  been  based  upon  and  took  its 
name  from  the  arithmetic-geometric  mean  inequality: 

T  T  §  T 

£  vt5t  i  rivt  1  if  vt,  6t  >  0  and  £  5t  =  L  (eqn  6.2) 

t=i  t=i  t=i 

The  equality  holds  only  when  v.  =  \'2=  If  ut  is  defined  as  ut=vt6t>  then 

Equation  6.2  becomes 

T  T  5 

Iut  >  n  (ut/6t)  1  (eqn  6.3) 

t=i  t=i 

The  equality  holds  if  =  u^/^Tuj..  Let  ut  be  a  posynomial  term  as  described  in 
Equation  6.4. 

ut  =  ctCPt^)]  =  ctn  xn&n’1  (eqn  6.4) 

A  posynomial  function,  f(x),  is  given  by  Equation  6.5. 

T 

fix)  =  S  ut  (ec*n  6-5) 

t=i 

When  the  posynomial  terms  are  substituted  into  Equation  6.3,  the  inequality  becomes 

T  T  n  a  6 

V  Uf  ^  n  {  Ccf  n  x  n,t]/6f}  1  (eqn  6.6) 

fci  1  t=i  1  n=x  n  V 

or 

I  ut  *  {FI  Cct/8t]6t}{n  xn<P} 

t=l  t=l  n=l 

where  <p  is  the  sum  over  t  of  an  t5t. 


(eqn  6.7) 


Since  the  only  restriction  on  St  has  been  that  they  be  positive  and  sum  to 
one,  they  may  be  chosen  such  that  (p  =  0  for  n  =  1,2, If  this  selection  is  made, 
Equation  6.5  becomes 


fix)  =  E  ctPtlx)  *  n  (ct/6t)  1 

t=i  t=i 

Since  equality  holds  when  5t  =  ut/Eut,  then 

t  T  6 

min  E  ut  =  max  n  (ct/5t)  1 


(eqn  6.8) 


(eqn  6.9) 


if  E  8t  =  1  and  Ean  t^t  =  0  for  n  =  1,2,...,N 

t=i  t=i  ’ 

Therefore,  the  minimization  of  the  posynomial  function  is  the  same  as  the 

maximization  of  the  nonlinear  function  in  Equation  6.9  subject  to  linear  constraints. 

*  * 

The  linearity  of  the  constraints  means  that  ot'  and  x  can  be  computed  with  linear 
algebra  as  explained  below,  instead  of  with  nonlinear  programming.  These 
minimization  and  maximization  problems  are  duals.  Since  equality  holds  if  and  only  if 
6t  =  ut/Eut,  it  follows  that  the  relationship  between  optimal  values  of  6  and  x  is 


St  =  {ctpt(X  )}/flX  )  or 


(eqn  6.10) 


*  *  f  * 
6,  -  {c,  n<x„  }  "-'ylfe  )• 


(eqn  6.1 1) 


If  the  degree  of  difficulty  is  zero,  then  the  matrix  of  exponents,  an  t,  with 

another  row  of  Is  appended  to  the  top  makes  a  square  matrix.  Rows  of  the  exponent 

matrix  correspond  to  variables  and  columns  correspond  to  terms.  The  row  of  l's 

corresponds  to  the  constraint  that  the  sum  over  t  of  5t  equals  1.  The  6  can  be 

* 

obtained  by  solving  a  set  of  T  simultaneous  linear  equations  Ao  =  b.  The  first 
element  of  the  b  vector  is  1  and  the  remaining  elements  are  0.  The  optimal  value  of 

♦  4: 

the  objective  function,  ffx  )  can  then  be  obtained  by  inserting  ot  into  Equation  6.12. 
Equation  6.12  is  based  upon  Equation  6.9. 


fix  )  =  n  (ct/«t )  1 

t=i 


(eqn  6.12) 


Finally,  the  optimal  values  of  the  decision  variables,  xn,  can  be  determined  by  solving 
the  set  of  T  equations  of  the  form  specified  in  Equation  6.13. 


E  an,tln(xn*)  =  ln[fl;x*)5t*/ct3  fort  =  1,2 . T. 


(eqn  6.13) 


This  set  of  equations  is  overconstrained  since  there  are  only  T-l  decision  variables. 

Therefore,  only  T-l  equations  are  required.  Solving  these  T-l  equations 

* 

simultaneously  produces  pn  =  )  which  are  then  converted  to  the  optimal  values 

*  Pn 

of  the  decision  variables  by  xn  =  e  . 
b.  Inequality  Constraints 

This  section  discusses  the  addition  of  posynomial  inequality  constraints  to 
the  unconstrained  problem  discussed  above.  For  notational  purposes  the  objective 
function  and  constraints  will  be  numbered  m=  0,1,2,3,...,M.  The  objective  function  is 
designated  m=0,  and  the  constraints  are  designated  m=  0,1,2, 3,...,M.  A  primal 
constrained  posynomial  would  have  the  form 


T  N 


Min  E  {c0,t  n  xn  °’n,t} 


t=l  n=l 


Subject  to: 


fraM  -  I  Vt  11  £  1  for  m  -  1,2 . M 


(eqn  6.14) 


(eqn  6.15) 


t=l  n=l 


where  xn  >  0,  n  =  1,2,3,... ,N.  If  6q  t  are  the  weights  for  the  terms  in  the  objective 
function,  then 


60,t  =  [c0,tP0,t^  )]/f0^  )  for  t=  1,2,3,. ..,T0. 


(eqn  6.16) 


If  X.m  are  the  Lagrange  multipliers  associated  with  constraint  m,  then 


I 

=  L  5m,t  and 

j.  —  i  ' 


(eqn  6.17) 


^m.t^rn  cm,t  Pm,t^^ 


(eqn  6. IS) 


(eqn  6.19) 


The  dual  geometric  program  is 


Subject  To: 


T 


t=i 

(eqn  6.20) 

M  T 

m 

I  I  am,n,t  5m,t  *  0  for  n  -  1,2,3 . N. 

m=I  t=l 

(eqn  6.21) 

T 

m 

“  I  «m,t 

t  =  l 

(eqn  6.22) 

and  6m,f  Xm  *  °- 

The  8m  t  are  calculated  using  Equations  6.20 

and  6.21  as  a  set 

simultaneous  linear  equations.  The  optimal  value  of  the  objective  function  is 


A-m 

calculated  by  multiplying  the  unconstrained  optimum  by  n(Xm)  111 

as  in  Equation 

6.23. 

r0(x*)  -  n  (ct/«t)6'  n  (X„/m 

t=l  m=l 

(eqn  6.23) 

* 

X  is  calculated  using  Equations  6.24  and  6.25. 

N  „  *  * 

Y  a0,n,tln^xn  )  =  W0(X  )50,t  /c0it)  for  1  =  U . TQ,  and 

n=l 

(eqn  6.24) 

N  *  *  * 

Y  am,n,t^n^xn  )  ~  ^n^m,t  '^cm,t^m  ) 

(eqn  6.25) 

for  t=  1,2,. ..,Tq;  m=  1,2, 3,... ,M;  and  Sm  t*>0. 

* 

As  with  the  unconstrained  problem  these  equations  are  linear  in  ln(xn  )  -  Pn  After 

solving  for  p  using  Equations  6.24  and  6.25  as  a  set  of  simultaneous  linear  equations, 

Pn 

the  decision  variables  are  calculated  using  xn  =  e  . 


C.  EXPLANATION  OF  VARIABLES 


•  B1(NT  +  1,NT*2)  holds  the  A  matrix  in  the  subroutine  which  solves  simultaneous 
linear  equations  of  the  form  Ax  =  b. 

•  B2(NT)  holds  the  b  vector  in  the  subroutine  which  solves  simultaneous  linear 
equations  of  the  form  Ax=  b. 

•  B3(NT,NT-1)  stores  the  exponents  of  the  variables  in  each  term.  Each  row  of 
the  matrix  corresponds  to  a  term;  each  column  corresponds  to  a  variable.  If  a 
variable  is  not  stated  explicitly  in  a  term,  then  its  entry  in  this  matrix  is  zero. 

•  CT(3,NC,MN)  holds  three  values  for  each  term.  CT(l,m,t)  holds  the  coefficient, 
cm  v  CT(2,m,t)  holds  the  weight,  Sm  t,  for  each  term.  CT(3,m,t)  holds  pm  t(x  ) 
for  each  term.  m=0  refers  to  terms  in  the  objective  function.  m=  1,2,...,NC 
refers  to  terms  in  the  mth  constraint.  t=  l,2,...,NT(m)  specifics  a  particular  term 
in  the  objective  function  or  a  constraint. 

•  FS  is  the  optimal  value  of  the  objective  function. 

•  11,12,  and  13  are  loop  counters. 

•  K2,K3,K4,...,K9  are  variables  in  the  simultaneous  linear  equation  solving 
subroutine.  This  subroutine  is  documented  in  Appendix  E. 

•  LM(NC)  holds  values  of  Xm,  m=  1,2,...,NC  where  Xm  is  the  sum  of  5m  t  for  the 
mth  constraint.  Xq  =  1. 

•  MN  is  the  maximum  number  of  terms  in  any  constraint  or  the  objective  function. 

•  NC  is  the  number  of  constraints. 

•  NT(NC)  is  the  number  of  terms  in  each  constraint. 

•  NT  is  the  number  of  terms  in  the  objective  function  and  the  constraints. 

•  NV  is  the  number  of  decision  variables,  i.e.  the  number  of  components  in  the 
vector  X- 

D.  INPUT 

Problem  parameters  are  entered  into  an  input  file,  GEOIN. DO,  before  the 
program  is  executed.  GEOIN. DO  must  contain  the  following  parameters  in  the  order 
specified. 

•  The  number  of  terms,  NT,  the  number  of  variables,  NV,  and  the  number  of 
constraints,  NC. 

•  The  number  of  terms  in  the  objective  function,  NT(0),  and  the  number  of  terms 
in  each  constraint,  NT(m)  m=  I,2,...,NC. 

•  The  coefficients  of  each  term,  cm  ,,  CT(l,m,t)  m  =  0,l,2,...,NC,  t=  l,2,...,NT(m). 


•  The  matrix  of  exponents,  an  t,  for  each  variable  in  each  term.  Row  n  of  the 
matrix  corresponds  to  the  variable  xn  n=  1,2,..., N.  Column  t  of  the  matrix 
corresponds  to  the  term  pt-(x),  t'=  1,2,. ..,T. 

The  input  file  for  the  problem  specified  in  Figure  6.1  is  in  Figure  6.2. 


Input  File: 

Probl em: 

4,  3,  1 

2  2 

Min  40x^2  +  2OX2X3 

At,  20 

Subject  To: 

.2,  .6 

1,0,  “1.  0 

1  1  -1 

.2xi"*X2‘^  +  . 6x2_^X3-2^3  ^  1 

0  1  0,  66666666667 

x1  >  0 

Figure  6.2  Sample  Input  File,  GEOIN. DO. 


E.  OUTPUT 

The  program  prints  the  following  output  to  the  screen. 

•  The  optimal  dual  variables,  $m  *. 

•  The  optimal  value  of  the  objective  function,  ffx  )• 

*\ 

•  The  value  of  each  pm  t(x  ' 

•  The  optimal  value  of  each  component  of  xn  . 

F.  EXPLANATION  OF  PROGRAM  COMPONENTS 

A  complete  program  listing  is  located  at  Appendix  D. 

1.  Initialization  And  Input,  Figure  6.3 

Line  110  opens  the  input  file,  GEOIN. DO.  Line  120  enters  the  number  of 
terms,  NT,  and  the  number  of  variables,  NV,  from  GEOIN. DO,  and  sets  K9  equal  to 
NT  for  use  in  the  simultaneous  linear  equation  solving  subroutine.  Line  122  checks 
whether  the  degree  of  difficulty  is  equal  to  zero.  If  it  is  not,  an  error  message  is  printed 
and  the  program  ends.  Line  130  enters  the  number  of  constraints,  NC,  and  dimensions 
the  vector  NT(NC),  which  holds  the  number  of  terms  in  the  objective  function  and  in 
each  constraint.  Line  140  enters  NT(m),  m=  0,1,2,. ..,NC,  and  computes  MN,  the 
maximum  over  m  of  NT(m).  Line  143  dimensions  the  matrices  required  for  the 
program.  Line  145  enters  the  coefficients  cm  t,  placing  them  in  CT(l,m,t). 

71 


100  'Gaonefric  Programming  Program 
110  OPEN"GEOIN"FORINPUTAS1 
120  INPUTP1 »NT  >NV : K9=NT 

122  IFNT-NV<>1THENPRINT"***ERR0R:  Dagreo  of  Difficulty  <>  0":END 
130  INPUTP1 >NC : DIMNT l NC ) 

140  MN=0 : FORIl=0TONC : INPUTP1 ,NT( 111 : IFNTI II  )>MNTHENMN=NT( II ) :NEXTI1 
143  0IMCT( 3 ,NC ,MN ) ,LM( NC ) ,B1< NT+1 ,NT*2  )  ,B2( NT  )  ,B3<  NT  ,NV ) 

14S  FORI1=OTONC:FORI2=1TONT<I1):INPUTS1,CT<1,I1,I2):NE>OT2:NEXTI1 
150  F0RI1=1T0NT1 0  ):B1< 1,11  )=1:NEXTI1 
155  F0RI1=NT<  0 )+lTONT:Bl(  1,11 )=0:NEXTI1 

160  F0RI1=2T0NT : F0RI2=1T0NT : INPUT#1 ,B11 11,12 ) :B3t 12 ,11-1 )=B1(  11,12 ) 
162  NEXTI2;NEXTI1 


Figure  6.3  Initialization  and  Input. 

Lines  150-162  enter  elements  of  the  A  matrix  required  to  solve  the 
simultaneous  linear  equations  A5  =  b  into  Bl(,).  Equations  6.20  and  6.21  are  the  basis 
for  this  set  of  simultaneous  linear  equations.  Lines  150-155  fill  the  first  row  of  B1  with 
ones  in  columns  corresponding  to  objective  function  terms  and  zeros  in  columns 
corresponding  to  other  terms.  Bl(l,t)  corresponds  to  Equation  6.20.  Lines  160-162 
enter  the  exponents  of  variables  in  each  term  into  B1  and  B3.  In  B1  rows  2  through 
NV  + 1  =  NT  correspond  to  variables  xn  and  columns  1  through  NT  correspond  to 
terms  t'  =  1,2,... ,T\  Storing  the  exponents  in  B3  is  necessary  because  the  simultaneous 
linear  equation  subroutine  changes  the  matrix  in  Bl,  and  the  exponents  are  required 
for  later  calculations.  B3  is  the  transpose  of  Bl  because  of  the  nature  cf  the 
calculation  for  which  B3  is  later  recalled. 

2.  Calculating  The  Weights  For  Each  Term,  Figure  6.4 


170  PRINT1"' : PRINT"**COMPUTING  DELTA'S**" 

172  B2( 1 )=1: FORIl=2TONT :B2( II )=0:NEXTI1 
180  GOSUB9800 

200  CLS: 11=1 : FORI2=OTONC : FORI3=lTONT( 12 ):CT( 2, 12,13 )=B1(  11,1) 

203  PRINT"DELTA< ">12>",">13l"  )  =  "j 

204  PRINTUSING"##### . ###*" >CT ( 2 ,12 ,13  ) : 11=11+1 : IFI1>5THENGOSUB600 

205  NEXTI3 :NEXTI2:GOSUB600:CLS 


m,r 


Figure  6.4  Calculating  5 


Line  172  places  the  simultaneous  linear  equation  b  vector  into  B2.  B2(l)=  1  is 
the  right  hand  side  of  the  constraint  in  Equation  6.20.  The  right  hand  sides  of  the 
constraints  based  upon  Equation  6.21  are  zero.  Line  180  calls  the  simultaneous  linear 
equation  solver  in  subroutine  9800  which  leaves  8  in  Bl(t',l),  t'=  1,2,. ..,T',  Lines 
200-205  place  8*  into  CT(2,m,t)  and  prints  8m  t  to  the  screen. 

3.  The  Optimal  Objective  Function  Value,  Figure  6.5 


210  PRINT . PRINT"**COMPLTTING  OPT  OBJ  FN  VALUE**" 

212  FORI1=OTONC : LM< II )=0: FORI2=lTONT< II ) :  LMl  II  )=LMt II  )+CT<  2,11 ,12  ) 
214  NEXTI 2 : NEXTI 1 
220  FS=1:FORI1=OTONC 

222  F0RI2=1T0NT( II ): FS=FS*(  CT(  1,11,12 )/CT(  2,11,12 ) )ACT( 2,11,12  ) 
224  NEXTI2:  FS=FS*(  LM<  II  )ALM(  II ) ): NEXTI  1 

229  PRINT""  :PRINT"F*  ="j :PRINTUSING"####.###" )FS:GOSUB600:CLS 


Figure  6.5  The  Optimal  Objective  Function  Value. 

Lines  212-214  compute  Xm  which  equal  the  sum  over  t  of  8mt  for  every 
constraint  and  the  objective  function.  Lines  220-224  compute  the  optimal  objective 
function  value  based  upon  Equation  6.23.  Line  229  prints  the  optimal  objective 
function  value. 

4.  Optimal  Decision  Variable  Values,  Figure  6.6 


230  'Compute  optimal  x(  n ) 

232  K9=K9-1 

234  FORIl=lTOK9: FORI2=lTOK9:Bl( 11,12 )=B3( 11,12  )  :NEXTI2:NEXTI1 

236  CC=1 : FORI1=OTONC : FORI2=lTONT( II ) 

237  CT(  3,11,12  )  =  <  CT<  2,11,12  )/CT<  1,11,12  >/LM<  II ) ) 

238  IFIl=0THENCTt 3,11,12  )=CT( 3,11,12  )*FS 

239  B2ICC  )=LOG( CT( 3 >11 ,12  )  ):CC=CC+1:NEXTI2:NEXTI1 

242  PRINT"P( m,t )*  =  opt.  value  of  term  t,  constr.  m,  divided  by  its  coef." 
244  FORI1=OTONC: FORI2=lTONT( II ) :PRINT"P(  "jll V" >")I2  >"  )*  ="> 

246  PRINTUSING"#### . ####” >CT<  3 , 1 1 ,1 2  ) : NEXTI2 : GOSUB600 : NEXTI 1 : CLS 
250  PRINT"": PRINT"**  Cotiputing  Opt  Values  Of  X(n)  **" 

260  GOSUB9800 

270  CLS : FORI 1=1T0K9:PRINT"X*( "jll j"  )  =  ”» 

272  PRINTUSING"######. #####" J EXP ( Bit  11,11) 

273  IFI1>5THENGOSUB600 
275  NEXTI1:GOSUB600:END 


Figure  6.6  Optimal  Decision  Variable  Values. 
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After  5  and  the  optimal  value  of  the  objective  function,  fQ(/  ),  have  been 


calculated,  this  section  calculates  x  •  The  basis  for  these  calculations  is  Equations  6.24 


and  6.25.  The  section  solves  a  set  of  simultaneous  linear  equations  A  [ln(x  )3  =  b  and 


then  solves  for  x  • 


Since  the  number  of  variables  is  one  less  than  the  number  of  terms,  line  232 
reduces  K9  by  one.  Line  234  creates  the  A  matrix  by  putting  the  matrix  of  exponents 
that  was  stored  in  B3  into  Bl. 

Lines  236-239  calculate  the  b  vector.  CC  is  a  counter  in  the  t'  subscripting 

system  which  controls  the  entry  of  b  vector  elements  into  B2(t').  Line  237  calculates 

the  portion  of  the  right  hand  side  that  is  common  to  all  terms.  Line  238  multiplies  that 

result  by  f()(x  )  for  objective  function  terms  which  produces  pt<(x  ).  Line  239  places 

ln[pt'(X*)J  into  B2(t').  Lines  242-246  print  Pm  t(X*)  t0  the  screen.  Line  260  calls  the 

’  * 

simultaneous  linear  equation  subroutine  which  solves  ACln(x  )]=b.  Lines  270-275 
print  x*  to  the  screen  and  end  the  program. 

5.  Subroutine  To  Stop  Screen  Printing,  Figure  6.7 


600  INPUT"**  Hit  ENTER  To  Continue:  "tZ9: RETURN 


Figure  6.7  Subroutine  To  Stop  Screen  Printing. 


Subroutine  600  is  used  to  interrupt  printing  loops  so  that  results  are  not 
scrolled  off  the  screen  before  the  operator  can  read  them. 

6.  Simultaneous  Linear  Equation  Subroutine,  Figure  6.8 

This  subroutine  solves  simultaneous  linear  equations  of  the  form  Ax=b.  K9 
is  the  dimension  of  the  square  A  matrix  and  the  x  and  b  vectors.  This  subroutine 
documented  in  Appendix  E,  The  Matrix  Algebra  Program. 


EXAMPLE  PROBLEMS 


1.  Example  #1 

A  design  engineer  wants  to  design  a  cylindrical  oil  storage  tank  with  a  storage 
capacity  of  lOOOtt  cubic  feet  to  put  on  an  existing  base.  If  the  cost  of  construction  is 


Sl/foot^  of  tank  surface,  what  are  the  optimal  dimensions  of  the  tank  and  how  much 
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'Simultaneous  Linear  Equation  Subroutine:  Ax=b 

'Invert  Matrix  A 

F0RK7=K9+1T02*K9: F0RK8=1T0K9 

IFK7=K8+K9THENB1( K8»K7  )=lELSEBlt K8 ,K7 1=0 

NEXTK8:N£XTK7 

F0RK7=1T0K9 

IFBK  K7 >K7  )*SGNI  Bl<  K7 »K7 )  )<1E -8THENGOSUB9910 

KZ-l/BK  K7,K7 ) :  F0RK6=1T02*K9:B1(  K7,K6  )=B1(  K7,K6  )*K2  :NEXTK6 

IFK7=K9THEN9865 

FORK8=K7+1TOK9:IFB1(K8,K7)=OTHEN9860 

K2=-B1IK8,K7) 

F0RK6=K7T02*K9:B1( K8.K6 )=B1( K8.K6 )+( K2*B1( K7.K6 ) ):NEXTK6 

NEXTK8 : NEXTK7 

F0RK7=K9T02STEP-1 

F0RK8=K7-1T01STEP-1 :  IFBK  K8  ,K7  )=0THEN9885 
K2=-B1(K8,K7) 

F0RK6=1T02*K9:B1IK8>K6  )=B1<  K8,K6  )+( K2*B1( K7»K6 ) ) :NEXTK6 

NEXTK8:NEXTK7 

'Mult  A  Inverse  by  b 

FORK7=lT0K9:BKK7,l)=0:F0RK8=lT0K9:Bl(K7,l)=Bl(K7,l)+BKK7,K8+K9)*B2(K8) 
NEXTK8 : NEXTK  7 : RETURN 
•Error  Routine 

IFERL>9700ANDERR=11THENPRINT"! ! IERROR:  Matrix  Is  Not  Invertable! ! !":ENO 
PRINT"Error  Code">ERR>''In  Line") ERL: END 
'SWITCH  ROWS 

F0RK5=K7+lT0K9:IFBKK5jK7  )#SGN(BKK5>K7 )  X1E-8THEN9940 
F0RK4=1T0K9*2 : K3=B1<  K7 ,K4 1  :B1<  K7 ,K4  )=B1<  K5 ,K4 ) 

B1 ( KS ,K4 ) =K3 : NEXTK4 : RETURN 

NEXTK5: PRINT "Error:  Matrix  Not  Invertable": END 


Figure  6.8  Simultaneous  Linear  Equation  Subroutine. 


will  it  cost?  The  tank  includes  the  cylindrical  siding  plus  the  top.  The  formulation  is  in 
Equations  6.26  and  6.27. 


Min  Sl(ttr^)  +  $l(27trh) 


(eqn  6.26) 


Subject  To: 


ttr^h  £  lOOOrt  r,h  >  0 


(eqn  6.27) 


The  objective  function  and  constraint  are  posynomials,  but  the  constraint  is  not  in  the 
<  1  form  required  by  the  program.  Putting  the  constraint  in  S  1  form  results  in 
Equation  6.28.  This  is  a  zero  degree  of  difficulty  problem  with  three  terms,  two 
variables,  one  constraint,  two  terms  in  the  objective  function,  and  one  term  in  the 


constraint. 


'ZSf’Xl 7 


>  V*  WV.WWV 


lOOOr’V1  £  1 


(eqn  6.28) 


The  coefficients,  ct',  t'  =  1,2,3  are  n,  2n,  and  1000  respectively.  The  input  file  and 
results  of  the  program  are  at  Figure  6.9. 


Input  File: 

Results: 

3,  2,  I 

2  1 

3.’ 142,  6.284,  1000 

0,  1,-1 

8l  =  1/3,  S2  =  2/3,  63  =  2/3 

f(X*)  =  $942.48 

2  1  -2 

Optimal  Values  for  r  =  10,  h  =  10. 

Figure  6.9  Input  File  and  Results  Of  Example  #1. 

The  economic  interpretation  of  8j  and  82  is  that  regardless  of  the  price  for 
steel,  1/3  of  the  cost  will  be  for  the  top  and  2/3  will  be  for  the  side.  If  the  top  and 
sides  were  constructed  of  different  types  of  steel  with  different  prices,  these  ratios 
would  not  change. 

2.  Example  #2 

This  problem  is  the  example  stated  in  Figure  6.1  The  input  file  is  in  Figure 
6.2.  8t-,  t'  =  1 ,2,3,4  are  .5,^5,  .5,  and  .75  respectively.  The  optimal  value  of  the 
objective  function  is  40.  The  optimal  values  of  xn,  n=  1,2,3  are  .5,  1,  and  1 
respectively. 


VII.  MATRIX  ALGEBRA  PROGRAM 


A.  GENERAL 

The  matrix  algebra  program,  MATALG,  performs  the  following  matrix  algebra 
functions:  matrix  addition,  multiplication,  and  inversion,  scalar  multiplication, 
calculation  of  determinants,  integer  exponentiation,  and  solutions  to  sets  of 
simultaneous  linear  equations.  MATALG  is  menu  driven.  The  main  menu  enables  the 
operator  to  enter  a  new  matrix,  print  the  answer  matrix,  or  call  one  of  the  functions 
listed  above.  Menus  produced  by  each  function  prompt  the  operator  for  required 
input.  Matrices  may  be  entered  from  the  M100  keyboard  or  from  a  RAM  file.  Output 
goes  to  the  M  100's  screen.  Intermediate  results  may  be  displayed.  Operations  are 
performed  in  the  conventional  left  to  right  order  in  which  matrix  operations  are  written 
out  on  paper,  e.g.  AxBxC*.  However,  if  the  series  of  operations  requires  altering 
that  order,  the  operator  may  store  one  matrix  for  future  recall.  Matrices  may  be 
entered  in  either  the  left  or  right  hand  position. 

B.  INPUT. 

The  matrix  input  subroutine  lets  the  operator  select: 

•  The  position9  of  the  matrix.  If  the  new  matrix  is  entered  for  the  left  side  of  the 
operation,  the  program  automatically  places  the  old  left  side  matrix  on  the  right 
side. 

•  Whether  the  matrix  will  be  entered  from  the  keyboard,  the  input  file 
MATIN. DO,  or  from  RAM  storage. 

•  Whether  the  matrix  will  be  scalar  multiplied  or  inverted. 

The  matrix  input  routine  must  be  accessed  from  the  main  menu  to  enter  the  first 
matrix.  From  then  on,  when  a  two  matrix  operation  is  selected  from  the  main  menu, 
the  subroutine  performing  that  operation  automatically  calls  the  input  subroutine  for 
the  second  matrix. 

If  the  input  file,  MATIN. DO,  is  used,  it  must  be  created  before  the  program  is 
run.  MATIN. DO  may  contain  more  than  one  matrix.  Matrices  must  be  preceded  in 
MATIN. DO  by  their  dimensions.  An  example  of  an  input  file  is  at  Figure  7.1. 

See  the  general  instructions  on  input  files  in  Chapter  2. 


9Lcft  or  right  side  of  the  operation. 
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Figure  7.1  Sample  Input  File,  MATIN. DO. 

When  matrices  are  entered  from  the  keyboard,  the  operator  will  be  prompted  for 
the  matrix  dimensions  and  for  each  matrix  element. 

C.  OUTPUT 

All  output  goes  to  the  screen  of  the  M100.  Output  may  have  up  to  three  digits 
to  the  left  and  four  digits  to  the  right  of  the  decimal  point.  If  this  configuration  is  not 
adequate,  the  operator  may  modify  the  format  at  the  line  numbers  specified  in  Figure 
7.2  for  the  corresponding  functions. 


LINE  FUNCTION 
1082  Determinant 

4075  Solution  to  Simultaneous  Linear  Equations 
6012  Other  Matrix  Output 

See  the  instructions  for  the  PRINT  USING  command  in  Reference  1. 


Figure  7.2  Line  Numbers  Of  Output  Formats. 

D.  DESCRIPTION  OF  VARIABLES 

•  A1(3,K,K)  holds  the  current  matrices.  A  1(1.,}  =  Left  side  matrix/primary 
matrix/current  intermediate  answer  matrix.  Al(2„)  =  Right  side/ secondary 
matrix.  Al(3„)  =  Matrix  being  stored. 

•  B1(K,K)  holds  the  answer  matrix  as  it  is  being  calculated. 

•  C(4)  and  R(4)  are  the  number  of  columns/rows  in  A1  or  Bl.  C{1)  and  R(l) 
correspond  to  A!{1„).  C(2)  and  R(2)  correspond  to  Al(2„).  C(3)  and  R(3) 
correspond  to  Al(3,,).  C(4)  and  R(4)  correspond  to  Bl. 

•  CC  is  the  row  counter  in  matrix  output  routine. 

•  CD  and  RD  are  the  column  and  row  of  element  to  be  changed. 


•  CH  is  the  selection  variable  for  the  main  menu. 

•  DET(2)  are  the  determinants  of  Al(l„)  and  Al(2„). 

•  EF  is  the  error  Flae.  0  -  No  terminal  error  has  been  made.  1  -  A  terminal  error 
has  been  made.  A  terminal  error  is  one  from  which  the  program  can  not  recover. 

•  FF  is  the  file  Flag.  1  -  MATIN'. DO  exists  in  RAM.  0  -  MATIN'. DO  does  not 
exist  in  RAM. 

•  II,  12,  13,  J 1 ,  J 2,  and  J3  are  loop  counters. 

•  K  is  the  largest  dimension  of  largest  matrix  to  be  processed. 

•  K4-K9  are  counters  in  the  simultaneous  linear  equation  subroutine. 

•  MF  is  the  multiplication/addition  flag.  0  -  Neither  the  multiplication  nor  the 
addition  subroutines  are  running.  1  -  The  multiplication  subroutine  is  running.  2 
-  Addition  subroutine  is  running. 

•  MI  is  the  matrix  indicator.  It  shows  which  matrix  in  A 1  ( 1  -3)  is  being  operated 
upon. 

•  MU  is  an  intermediate  multiplier  in  the  determinant  and  matrix  inversion 
subroutines. 

•  OF  is  the  output  flag.  1  -  Send  output  to  screen.  0  -  Suppress  output  to  screen. 

•  SF  is  the  simultaneous  linear  equation  (SLE)  flag:  1  -  The  SLE  subroutine  is 
running.  0  -  The  SLE  subroutine  is  not  running. 

•  XP  is  the  umber  of  times  the  matrix  will  be  multiplied  times  itself  in  the  integer 
exponentiation  subroutine. 

•  Z9  is  a  general  purpose  variable. 


E.  DESCRIPTION  OF  PROGRAM  COMPONENTS 
A  complete  program  listing  is  located  at  Appendix  E. 
1.  Initialization,  Figure  7.3 


100  CLS: PRINT"": PRINT"  ***  MATRIX  ALGEBRA  PROGRAM  ***": PRINT"" 

105  PRINT"IS  INPUT  MATRIX,  ’MATIN. 00'  IN  RAM?" : INPUT"  0=NO>  1=YES"jFF 

107  IFFF=1THEN0PEN"MATIN"F0RINPUTAS1 

300  PRINT"*#Enter  The  Single  Largest  Dimension  of" 

305  INPUT"The  Largest  Matrix  To  Be  Processed:  " »K 

310  DIMAll 3 ,K,K  )  >B1<  K+l ,K*2  )  ,R<  4 ) ,C<  4  )  >DET( Z  ) : MI  =  1 :OF=l :SF=0 


Figure  7.3  Initialization  Section. 

Line  100  prints  the  program  title.  Line  105  permits  the  operator  to  specify 
whether  file  MATIN. DO  will  be  used  as  a  source  of  input  and  sets  the  file  flag,  FF, 
accordingly.  Line  107  opens  MATIN. DO  for  input  if  it  is  to  be  used. 
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Lines  300-305  require  the  operator  to  specify  the  largest  dimension,  K,  of  the 
largest  matrix  to  be  processed.  Line  310  dimensions  matrices  A1  and  B1  and  vectors 
R,  C,  and  DET  and  initializes  flags  MI,  OF,  and  SF. 

2.  Main  Menu,  Figure  7.4 

This  section  permits  the  operator  to  select  the  next  major  operation  to  be 
conducted  and  calls  the  subroutine  performing  that  operation. 


501  CLS:EF=0:PRINT"****MATRIX  ALGEBRA  PROGRAM  MENU****" 

504  PRINT"  1.  Enter  Starting  Left  Side  Matrix” 

505  PRINT"  2.  Matrix  Inversion" 

506  PRINT"  3.  Matrix  Addition" : PRINT"  4.  Matrix  Multiplication" 

508  PRINT"  5.  Simultaneous  Linear  Equations" 

509  PRINT"  6.  Print  Current  Answer  Matrix" : PRINT"  7.  Other  Options" 

510  INPUT"  **Enter  Number:  "sCH 

512  IFCH=1THENMI=1 :Z9=0 : GOSUB7006 

513  IFCH=2THENMI=1:GOSUB2000 

514  IFCH=3THENGOSUB3000 

515  IFCH=4THENGOSUB5000 

516  IFCH=5THENGOSUB4000 

517  IFCH=6THENGOSUB6000 

518  IFCHO7THENG0T0501 

520  CLS: PRINT"***MORE  CHOICES PRINT"  1.  Determinant" 

524  PRINT"  2.  Matrix  Integer  Exponentiation" 

526  PRINT"  3.  Store  Current  Matrix" 

530  PRINT"  4.  Retrieve  Stored  Matrix" : PRINT"  5.  Scalar  Multiplication 
532  PRINT"  6.  Other  Options" : INPUT"**Enter  Number:  "»  CH 
540  IFCH=1THENMI=1:GOSUB1000 

548  IFCH=2THENMI=1:GOSUB7600 

549  IFCH=3THENG0SUB8000 

550  I FCH=4THENMI=1 : GOSUB8Z00 
560  IFCH=5THENMI=1:GOSUB5100 
570  GPT0501 


Figure  7.4  Main  Menu. 

Line  501  initializes  EF  and  prints  the  main  menu  header  to  the  screen.  Lines 
502-510  print  the  First  screen  of  options  and  prompt  the  operator  for  a  selection.  Lines 
512-518  call  the  subroutine  selected  by  the  operator  or  branch  to  the  second  screen  of 
options.  Lines  520-532  print  the  second  screen  of  options  and  prompt  the  operator  for 
a  selection.  Lines  540-570  call  the  subroutine  selected  by  the  operator  or  return  to  the 
first  screen  of  options. 

3.  Pause  Control  Subroutine,  Figure  7.5 

Lines  700-702  stop  the  program  to  permit  the  operator  to  view  material  on  the 
screen  and  permit  continuation  by  pressing  the  ENTER  button. 


A*"  V,  -V*^,  *\  *  .  *%  v\.  a/.  /,  v. 


700  'PAUSE  CONTROL 

702  INPUT"**  Hit  ENTER  To  Continue" >Z9: RETURN 


Figure  7.5  Pause  Control  Subroutine. 
4.  Modifying  the  Secondary  Matrix,  Figure'7.6 


800  'INTERMEDIATE  MODIFICATIONS 
810  PRINT"**Modify  The  2nd  Matrix?" 

812  INPUT"  0=No,  l=Invert,  2=Scalar  Multiply 
815  IFZ9=0THENRETURN 

820  MI=2:IFZ9=1THENGOSUB2000ELSEGC)SUB5100 
825  GOT0810 


Figure  7.6  Modifying  The  Secondary  Input  Matrix. 

This  subroutine  permits  the  operator  to  invert  the  second  matrix  of  a  two 
matrix  operation  or  multiply  that  matrix  by  a  scalar.  Lines  810-812  print  the  options 
to  the  screen  and  prompt  the  operator  for  a  selection.  Line  815  causes  a  return 
without  the  matrix  being  modified  if  appropriate.  Line  820  sets  the  matrix  indicator, 
MI,  to  two  and  calls  the  matrix  inversion  or  scalar  multiplication  subroutine.  Line  825 
starts  the  subroutine  again,  permitting  the  operator  to  select  another  option. 

5.  Determinant  Calculation,  Figure  7.7 

Lines  1005-1008  test  whether  the  matrix  is  square  and  print  an  error  message 
if  it  is  not.  Line  1010  copies  the  matrix  to  be  inverted  into  B1  where  the  calculations 
will  be  conducted.  Line  1020  initializes  the  value  of  the  determinant  as  one  and  the 
row  counter,  II. 

Line  1021  checks  whether  the  diagonal  element  in  the  current  row,  Rc,  is  zero. 
If  so,  then  the  row  switching  subroutine  is  called.  If  all  the  rows  below  Rc  have  0's  in 
column  II  then  the  determinant  is  zero.  If  a  non-zero  element  can  be  found  below  Rc 
in  column  II  then  that  row  is  switched  with  Rc  and  the  determinant  is  multiplied  by  -1. 
Line  1022  tests  for  a  termmal  error  from  the  row  switching  subroutine.  Line  1023 
multiplies  the  determinant  by  diagonal  element  II  and  branches  to  the  end  of  the 


1000  'CALC  DETERMINANT 
1005  IFR( MI  )=C( MI JTHEN1010 

1007  PRINT"ERROR:  Number  of  rows/columns  not  equal 

1008  PRINT”  MATRIX  IS  NOT  INVERTABLE  !*' :  G0SLIB 700:  EF=1 :  RETURN 

1010  F0RI1=1T0R<  MI ) : F0RI2=1T0C( MI  )  :B11 II ,12  )=A1( MI ,11 ,12  1  :NEXTI2 :NEXTI1 

1020  0ETIMI  )=1:F0RI1=1T0R(MI ) 

1021  IFB1(  II  ,11  XSGN1  Bl(  11,11 )  X1E-10THENGOSUB1900ELSE1023 

1022  IFEF=1THEN1008 

1023  DET( MI )=DET<  MI )*B1<  11,11 ) :IFI1=R( MI 1THEN1080 

1025  F0RI3=1T0C(MI )  :B1(  11,13  )=B1(  1.1 ,13  )/Bll  11,11 )  :NEXTI3 
1030  F0RI2=I1+1T0R(MI ) : IFB1( 12,11 X0THEN1060 

1040  F0RI3=I1T0C(MI  ):B1(  12,13  )=B1(  12, 13  MBit  12,11  )*B1<  11,13  )  >:NEXTI3 
1060  NEXTI2:NEXTI1 

1080  IFOFOITHENRETURN 

1081  PRINT“#*Det.  Of  Matrix  "}MI>”  Is:  "j 

1082  PRINTUSING"##### .  #####"  >  DET  ( MI ) :  GOSUB700 
1090  RETURN 


Figure  7.7  Determinant  Calculation. 

subroutine  if  the  Rc  is  the  last  row.  Line  1025  divides  all  elements  in  Rc  by  the 
diagonal  element  in  Rc>  Lines  1030-1060  update  the  elements  in  Rc  and  below  in  in 
column  II  and  to  the  left.  Line  1080  tests  whether  output  is  to  be  printed  to  the 
screen.  -If  not,  the  subroutine  ends.  If  so,  lines  1081-1082  print  the  determinant. 

6.  Row  Switching  Subroutine,  Figure  7.8 


1900  'SWITCH  ROWS 

1910  FORJ=Il+lTOR<  MI ) :  IFB1(  J  ,11  )*SGNf  Bit  J  ,11 )  X1E-10THEN1940 
1920  FORJl=lTOC(  MI  )*2:TE=B1<  II, J1 )  :B1(  II, J1  XBll  J,J1 ) 

1930  Bl( J,J1XTE:NEXTJ1:GOT01950 

1940  NEXTJ:EF=0: RETURN 

1950  DETIMI  )=-DET( MI  ): RETURN 


Figure  7.8  Row  Switching  Subroutine. 

This  subroutine  is  called  when  the  determinant  or  inversion  subroutines  try  to 
pivot  on  a  row,  Rjj,  with  a  zero  in  the  main  diagonal  element.  The  subroutine  looks 
for  the  first  row  below  Rj|  that  has  a  nonzero  element  in  column  II.10  Line  1910 
searches  the  rows  below  Rjj  for  a  row  with  a  nonzero  element  in  the  appropriate 


10The  same  column  as  the  zero  on  the  main  diagonal  in  Rj  j. 


column.11  Lines  1920-1930  switch  the  elements  of  the  rows  using  TE  as  an  intermediate 
storage  variable.  If  none  of  the  rows  below  Rjj  have  a  nonzero  element  in  the 
appropriate  column,  then  the  matrix  is  not  invertable  and  EF  is  set  to  one  in  line  1940. 
If  a  row  switch  was  made,  the  determinant  changes  sign  in  line  1950. 

7.  Matrix  Inversion  Subroutine,  Figure  7.9 


2000  'MATRIX  INVERSION 

2010  OF=0 : GOSUBIOOO : I FDET ( MI )*SGN(DET(MI ) )>1E-100REF=1THEN2017 

2015  PRINT"*ERROR:  Determinants .  MATRIX  NOT  INVERTABLE !":GOSUB700:EF=1 

2017  I FE  F=1THENRETURN 

2020  F0RI1=1T0R(MI ) : F0RI2=1T0C(MI ):B1( 11,12 )=A1(MI ,11,12 ) :NEXTI2:NEXTI1 
2030  F0RI1=C< MI )+lT02*C(MI >: F0RI2=1T0R< MI ) 

2032  IFI1=I2+R( MI  )THENB1( 12,11  )=1ELSEB1( 12,11  )=0 
2035  NEXTI2:NEXTI1 
2040  F0RI1=1T0C( MI ) 

2045  IFB1(  II  ,11  )*SGN(  61(11,11)  X1E-10THENGOSUB1900ELSE2055 

2046  IFEF=1THEN2015 

2055  MU=1/B1( 11,11): FORI3=1T02*C( MI ) :B1( 11,13  )=B1( II ,13  )*MU:NEXTI3 
2057  IFI1=C(MI JTHEN2080 

2060  FORI 2=1 1+1TORI MI )  :IFB1( 12,11 )=0THEN2075 
2065  MU=-B1( 12,11) 

2070  F0RI3=I1T02*C(MI  ):B1( 12,13  )=B1( 12,13 )+( MU*B1( 11,13 )  ):NEXTI3 

2075  NEXTI2:NEXTI1 

2080  F0RI1=C( MI  )T02STEP-1 

2100  F0RI2=I1-1T01STEP-1:IFB1( 12,11  )=0THEN2130 
2110  MU=-B1( 12,11) 

2120  F0RI3=1T02*C(MI ):B1( 12,13 )=B1( 12,13 )+( MU*Bll 11,13 ) ) :NEXTI3 

2130  NEXTI2 :NEXTI1 

2140  F0RI1=1T0C( MI ) : F0RI2=1T0R( MI ) 

2145  All  MI ,12,11 )=B1< 12 ,11+C( MI ) ) :NEXTI2:NEXTI1 
2190  OF=l: RETURN 


Figure  7.9  Matrix  Inversion  Subroutine. 

The  matrix  inversion  subroutine  places  the  matrix  to  be  inverted,  p,  and  an 
identity  matrix,  I,  into  BI.  Each  row  of  B1  holds  a  row  of  p  and  the  corresponding 
row  from  I.  Elementary  row  operations  are  conducted  on  p  in  B I  to  change  that 
portion  of  Bl  to  an  identity  matrix.  The  same  elementary  row  operations  are 
conducted  on  the  portion  of  Bl  that  started  as  an  identity  matrix.  When  the  portion 
of  Bl  which  started  as  p  is  changed  to  I,  then  the  portion  of  Bl  which  started  as  I 
becomes  p'1. 


^The  decision  rule  actually  looks  for  an  element  outside  the  range  0  ±  10'^. 


Line  2010  stops  intermediate  results  from  being  printed  to  the  screen  by 
setting  OF  to  zero.  Line  2010  also  calls  the  determinant  calculation  subroutine  and 
tests  whether  the  determinant  is  equal  to  zero.  If  the  determinant  equals  zero,  then 
lines  2015-2017  print  an  error  message,  set  the  error  flag,  and  terminate  the  inversion 
subroutine.  Line  2020  copies  p  into  Bl.  Lines  2030-2035  place  an  identity  matrix  with 
the  same  dimensions  as  p  into  B 1  with  p. 

Lines  2040-2075  Line  2045  checks  whether  the  pj  j  j  j  is  zero.  If  so,  then  the 
row  switching  subroutine  is  called.  If  the  row  switching  subroutine  can  not  find  a  row 
for  which  element  II  does  not  equal  zero,  then  the  matrix  is  not  invertable  and  line 
2046  branches  to  the  error  message.  Line  2055  calculates  the  constant,  MU,12  and 
multiplies  row  II  by  MU.  Since  lines  2060-2075  do  not  apply  to  the  last  row  of  p,  line 
2057  branches  around  them  if  1 1  points  to  the  last  row. 

Lines  2060-2075  perform  the  elementary  row  operations  to  change  to  zero  the 
elements  of  column  II  that  are  in  rows  below  row  il.  Lines  2080-2130  perform  the 
elementary  row  operations  which  change  to  zero  the  elements  of  p  above  the  main 
diagonal.  Lines  2140-2145  copy  the  inverted  matrix  from  Bl  back  to  the  appropriate 
section  of  Al.  Line  2190  turns  the  output  back  on  by  resetting  OF  and  terminates  the 
subroutine  with  a  return. 

8.  Matrix  Addition,  Figure  7.10 

3000  'MATRIX  ADDITION 

3010  MF=2 :  GOSCJB7000 :  GOSUBSOO :  IFEF=1TH£NRETURN 

3015  FORIl=lTOR(  1 ) ;  FORI2  =  lTOC(  1):  All  1,11,12  )=A1(  1,11,12  )+Al<  2,11,12) 

3020  NEXTIZ : NEXTI 1 : GOSUB6000 : MF =0 : RETURN 

Figure  7.10  Matrix  Addition  Subroutine. 

Line  3010: 

•  Sets  MI  =  2  indicating  to  the  input  subroutine  that  it  is  being  called  from  the 
matrix  addition  subroutine. 

•  Calls  subroutine  7000  to  enter  the  second  matrix. 


^Multiplying  row  II  by  MU  makes  main  diagonal  element  Pii  ti  equal  to  one, 
i.e.  its  identity  matrix  value.  ’ 


•  Calls  subroutine  800  to  permit  the  operator  to  invert  the  second  matrix  or 
multiply  it  by  a  scalar. 

•  Evaluates  whether  a  terminal  error  was  made  in  either  subroutine  7000  or  800 
and,  if  so,  terminates  the  matrix  addition  subroutine. 

Lines  3015-3020  add  the  elements  of  Al(l„)  and  A2(2„).  Line  3020  also  calls 

subroutine  6000,  printing  the  answer,  resets  MF  =  0,  and  terminates  the  matrix  addition 

subroutine. 

9.  Simultaneous  Linear  Equations,  Figure  7.11 


4000  'SIMULTANEOUS  LINEAR  EQUATIONS 

4010  CLS : PRINT"**Solves  Ax=b.  Choices PRINT"  I.  Enter  b  Vector" 

4012  PRINT"  2.  Change  An  Element  In  Matrix  A" 

4013  PRINT"  3.  Solve  Current  Ax=b" 

4014  PRINT"  4.  Return": INPUT"  *  Select  A  Number:  "jCC 
4020  IFCC=1GOT04040 

4022  IFCC=2GOT04050 
4024  IFCC=3GOTO4060 
4026  IFCC=4THEN  RETURN 
4035  GOTO  4000 

4040  MI=2:R( 2 )=C< 1 ):Cl 2 )=l:GOSUB7040:G0T04000 

4050  INPUT"**Row,  Column  Of  Matrix  A  To  Be  Changed:  ">RD,CD 

4052  PRINT"  -  Enter  Row"jRDl",  Column" iCD } " : " J : INPUTA1 ( 1 , RD  >CD  ) : GOT04000 

4060  MI =1 : SF=1 : 0F=0 : GOSU88000 : GOSUB20OO 

4064  IFEF=0THEN4070 

4065  PRINT"*Solution  Not  Uniquely  Determinable" :GOSUB700: RETURN 
4070  GOSU85000:CC=0:FORIl=lTOR(2):CC=CC+l:PRINT"x("}Il»")  =  "» 

4075  PRINT  USING"##### .  ####"  lBl(  II ,  1 ) :  IFCO6THENG0SUB700 :  CC=0 
4080  NEXTI 1 : GOSUB700 : SF=0 : GOSUB8200 : GOT04000 


Figure  7.11  Simultaneous  Linear  Equation  Solving  Subroutine. 


This  subroutine  solves  sets  of  linear  equations  of  the  form  Ax=b  where  A  is 
m  by  m  matrix  of  rank  m  and  x  and  b  are  vectors  of  length  m.  Matrix  A  must  be 
entered  as  the  primary  matrix  before  this  subroutine  is  called.  The  subroutine  prompts 
the  operator  to  enter  b.  The  subroutine  also  permits  the  operator  to  change  individual 
elements  of  the  A  matrix. 

Lines  4010-4014  print  a  header  and  a  menu  of  options  and  prompt  the 
operator  to  select  an  option.  Lines  4020-4035  transfer  control  to  execute  the  option 
selected.  Line  4040  sets  MI  =  2,  indicating  that  b  will  be  stored  in  Al(2„),  dimensions 
the  b  vector,  and  calls  subroutine  7040  to  input  b.  Line  4050  prompts  the  operator  for 
the  row  and  column  of  the  element  in  A  to  be  changed.  Line  4052  prompts  the 
operator  to  enter  the  new  value  for  that  element. 


Lines  4060-4080  solve  the  system  of  equations.  Line  4060  sets  MI,  SF,  and 
OF  and  calls  subroutines  which  store,  then  invert,  the  A  matrix.  Lines  4064-4065  print 
an  error  message  if  the  A  matrix  is  not  invertable.  Lines  4070-480: 

•  Multiply  A**  by  b  producing  the  solution  vector,  x. 

•  Print  the  solution  vector. 

•  Reset  SF  =  0. 

•  Retrieve  the  stored  A  matrix. 

10.  Matrix  Multiplication  Subroutine,  Figure  7.12 

5000  'MATRIX  MULT 

5010  MF=1:IF  SF=1THEN5020 

5015  MI=2:GOSU87000:GOSUB800: IFEF=1THENRETURN 

5020  R( 4 )=R< 1 )  :C( 4  )=C<  2  ) : F0RI1=1T0R<  4  ) : FORI2=lTOC( 4 ) :B1< 11,12 )=0 
5022  FORI3=lTOCl  1 ) :B1( II ,12 )=A1< 1,11,13  )*A1( 2,13,12  J+Bll 11,12  ) 

5024  NEXTI3:NEXTI2:NEXTI1:MF=0 
5050  IF  SF=OTHENGOSU87SOO : GOSUB6000 
5060  RETURN 

Figure  7.12  Matrix  Multiplication  Subroutine. 

Line  5010  sets  MF13  and  branches  to  avoid  the  input  subroutines  in  line 
5015  if  SF  equals  one.14  Line  5015  calls  the  subroutines  to  enter  and  modify  the  second 
matrix.  Lines  5020-5024  set  the  dimensions  B1  and  perform  the  multiplications  and 
additions  required  to  place  the  answer  matrix  in  Bl.  If  SF  equals  zero,15  then  line 
5050  calls  subroutines  which  copy  the  answer  from  Bl  to  Al(l„)  and  print  the  answer 
matrix. 

1 1.  Scalar  Multiplication  Subroutine,  Figure  7.13 

This  subroutine  multiplies  matrix  A1(MI„)  by  a  scalar,  SM.  Lines 
5110-51 15  prompt  the  operator  to  enter  the  scalar  and  conduct  the  multiplication. 


13MF=  1  indicates  that  matrix  multiplication  is  being  performed. 

14That  is,  if  the  matrix  multiplication  subroutine  is  called  from  the  simultaneous 
linear  equation  subroutine. 

15That  is,  this  subroutine  is  not  being  called  from  the  simultaneous  linear 
equation  subroutine. 


S100  'SCALER  MULT 

5110  INPUT"Enter  Scalar  Multiplier :"jSM: FORIl=lTORt MI ) : F0RI2=1T0C< MI ) 
5115  Al<  MI >11,12 )=A1(MI>I1>I2 )*SM:NEXTI 2 :NEXTI1: RETURN 


Figure  7.13  Scalar  Multiplication  Subroutine. 
12.  Subroutine  To  Print  Al(l„)>  Figure  7.14 


6000  'PRINT  OUTPUT  MATRIX 

6010  PRINT"  **  Current  Answer  Matrix:*':CC=0: FORI1=1TOR(  1 ) :CC=CC+1 
6012  F0RI2=1T0C( 1 ): PRINTUSING"####.####" »A1( 1>I1 ,12 )  > :NEXTI2:PRINT“" 
6050  IFCC=3THEN60SUB700:CC=0 
6070  NEXTI1:60SUB700: RETURN 


Figure  7.14  Subroutine  To  Print  The  Primary  Matrix. 

Lines  6010-6070  print  a  header  and  then  print  Al(l„).  The  matrix  is  printed 
up  to  three  rows  at  a  time. 

13.  IVIatrix  Input  Subroutine,  Lines  7000-7050 
a.  Input  Matrix  Configuration,  Figure  7.15 

The  section  in  Figure  7.15  prompts  the  operator  to  specify: 

•  Whether  the  matrix  to  be  entered  will  go  on  the  left  or  right  hand  side  of  the 
operation. 

•  Whether  the  matrix  will  be  entered  from  the  keyboard,  MATIN. DO,  or  retrieved 
from  RAM  storage.  The  dimensions  of  the  incoming  matrix  are  entered  from  the 
appropriate  source. 

Line  7001  prompts  the  operator  to  specify  whether  the  incoming  matrix  will 
be  on  the  left  or  right  side  of  the  operation.  If  the  incoming  matrix  is  to  be  on  the  left 
side,  lines  7003-7004  move  the  matrix  in  A I(  1  „)  to  Al(2„).  Lines  7006-7008  print  an 
appropriate  header.  Lines  7009-7011  print  the  source  options  for  the  incoming  matrix. 
The  option  to  enter  a  matrix  from  MATIN. DO  will  be  printed  only  if  MATIN. DO  has 
been  created  (FF  =  1)  and  the  end  of  MATIN. DO  has  not  been  reached  (EOF(1)  =  0). 
Lines  7012-7013  transfer  control  to  enter  the  matrix  from  the  appropriate  source. 
Lines  7014-7018  enter  the  dimensions  of  the  new  matrix  from  the  appropriate  source. 


7000  ‘MATRIX  INPUT 

7001  CLS : PRINT"" : PRINT "Will  Thi*  Matrix  Be  0n:“:INPUT"  0=Left,  l=Right"jZ9 

7002  IFZ9=1THEN7006 

7003  R(  2  )=R(  1 )  :C(  2  )=C(  1 ): F0RI1=1T0R( 1 ) : F0RI2=1T0C( 1):A1(2,I1,I2 )=A1( 1,11,12 ) 

7004  NEXTI2:NEXTI1:MI=1 

7006  CLS:IFMI=2THEN7008 

7007  PRINT"«*Choices  For  Left  Hand  Matrix" :GOT07009 

7008  PRINT"**Choices  For  Right  Hand  Matrix:" 

7009  PRINT"  1.  Enter  Matrix  From  Keyboard” 

7010  PRINT"  2.  Retrieve  Stored  Matrix" :IFFF<>1THEN7012 

7011  IFE0F( 1 )=0THENPRINT"  3.  Enter  Matrix  From  MATIN. DO" 

7012  INPUT“*#Enter  A  Number:  “»Z9:IFZ9=1THEN7015 

7013  IFZ9=3THEN7018 

7014  R(MI  )=R( 3 ) :C( MI  )=CI 3 ) :GOTO7020 

7015  PRINT”  **Enter  The  Rows,  Columns" 

7017  INPUT“In  The  Next  Matrix:  “jRCMI  ),C(MI J:GOT07020 

7018  INPUT«1,R(MI),C(MI) 


Figure  7.15  Input  Matrix  Configuration. 


b.  Detection  Of  Dimensioning  Errors ,  Figure  7.16 


7020  IF  MFO1THEN7030 

7021  IFR( 2  )=C< 1 1THEN7030 

7022  PRINT"**ERROR :  Columns  in  LEFT  MATRIX  ="jC( 1) 

7024  PRINT"  Rows  In  Right  Matrix  =">R(2) 

7026  PRINT"Theso  Must  Be  Equal  For  Matrix  Mult? f":GOSUB700:EF»l:GOT07006 

7030  IF  MFO2THEN7035 

7031  IF(  R(  1  )=R(  2  )ANDC(  1  )=C(  2  )  JTHEN  7035 

7032  PRINT"**ERROR :  Diinens  ions  For-Both  Input" 

7034  PRINT"Ma trices  Must  Be  Equal!!":GOSUB700:EF=l:GOT07006 


Figure  7.16  Detection  Of  Dimensioning  Errors. 

If  the  matrix  being  entered  is  the  second  matrix  in  a  matrix  multiplication 
operation,  then  lines  7020-7026  check  whether  the  number  of  columns  in  the  left  matrix 
is  equal  to  the  number  of  rows  in  the  right  matrix.  If  not,  then  an  error  message  is 
printed  and  control  is  transfered  to  the  beginning  of  the  matrix  input  subroutine.  If 
the  matrix  being  entered  is  the  second  matrix  in  a  matrix  ar  dition  operation,  then  lines 
7030-7034  check  whether  the  dimensions  of  the  left  and  righ .  matrices  are  the  same.  If 
not,  then  an  error  message  is  printed  and  control  is  transfe  ed  to  the  beginning  of  the 
matrix  input  subroutine. 
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c.  Matrix  Input  Section ,  Figure  7.17 


7035  I f Z9=2THENGOSUB8200 : RETURN 
7035  IFZ9=3THEN7050 

7037  PRINT"  **Fill  Matrix  Row  By  Row:": PRINT"" 

7040  F0RI1=1T0R( MI ) : F0RI2=1T0C( MI ) 

7042  PRINT" -Enter  Row">Il»"And  Column"a2»":"j 
7044  INPUTA1 ( MI , I 1 , 12 ) : NEXTI2 : PRINT "" : NEXTI 1 : RETURN 

7050  FORI 1= 1TOR l MI ) : FORI 2 =1T0C ( MI ) : INPUT#1 , A1 ( MI , I 1 , I 2 ) : NEXTI 2 : NEXTI 1 : RETURN 


Figure  7.17  Matrix  Input  Section. 

If  the  incoming  matrix  is  to  be  retrieved  from  RAM  storage,  line  7035  calls 
the  appropriate  subroutine.  Line  7037-7044  enter  the  incoming  matrix  from  the 
keyboard.  If  the  incoming  matrix  is  to  be  entered  from  MATIN. DO  then  line  7036 
transfers  control  to  7050  which  performs  the  entry. 

14.  Copy  B1  Into  A  1(1„),  Figure  7.18 


7500  'COPY  B1  INTO  Al(l,,) 

7510  R( 1 )=R( 4 ) :C<  1 )=C( 4 ) : F0RII=1T0R( I ) : F0RI2=1T0CI I ) 
7512  AK1, II, 12  )=B1<  11,12  ): NEXTI 2: NEXTII:  RETURN 


Figure  7.18  Subroutine  to  copy  B1  into  Al(l„). 

Lines  7510-7512  dimension  al(l„)  and  copy  B1  into  Al(l„). 
15.  Matrix  Integer  Exponentiation,  Figure  7.19 


7600  'MATRIX  INTEGER  EXPONENTIATION 

7610  CLS: PRINT""  :INPUT"**Enter  Integer  Exponent  >  2:  ")XP 
7620  R( 2  >=RI 1 1  :C( 2  )=CI 1 1 : F0RI1=1T0R( 1 ) : F0RI2=1T0CI 2  ) 

7622  Al<  2,11,12 )=A1(  1,11,12 ):NEXTI2:NEXTI1 

7630  SF=1 : FOREX=2TOXP : GOSUB5020 : GOSUB7500 : NEXTEX : GOSl»6000 : SF=0 : RETURN 


Figure  7.19  Matrix  Integer  Exponentiation  Subroutine. 


This  subroutine  raises  the  primary  matrix  to  an  integer  power  greater  than  or 
equal  to  two.  Line  7610  prompts  the  operator  to  enter  an  exponent,  XP.  Lines 
7620-7622  dimension  Al(2„)  and  copy  Al(l„)  into  A2(2„).  Line  7630: 

•  Sets  SF=  1.  This  suppresses  printing  of  intermediate  results. 

•  Performs  XP-1  matrix  multiplications. 

•  Prints  the  final  result. 

16.  Storage  and  Retrieval  Subroutines,  Figure  7.20 
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Figure  7.20  Storage  and  Retrieval  Subroutines. 

Lines  8010-8012  dimension  Al(3„)  and  store  A  1(1  „)  in  Al(3„).  Lines 
8210-8212  dimension  Al(l„)  and  retrieve  Al(l„)  from  Al(3„). 

F.  SIMULTANEOUS  LINEAR  EQUATION  SUBROUTINE,  FIGURE  7.21 

Many  programs  require  a  simultaneous  linear  equation  solver.  Often  these 
programs  compute  the  A  matrix  as  part  of  the  program  and  use  the  results  in 
subsequent  calculations.  The  following  subroutine  may  be  inserted  in  other  programs 
without  requiring  the  loading  of  the  entire  matrix  algebra  program. 

This  subroutine  follows  the  same  algorithm  as  the  simultaneous  linear  equation 
subroutine  in  the  Matrix  Algebra  Program.  K9  is  the  dimension  of  the  A  matrix.  B1 
holds  the  A  matrix;  B2  holds  the  b  vector.  The  x  vector  is  returned  in  the  first  column 
of  Bl.  If  the  A  matrix  is  to  be  used  later,  it  must  be  stored  somewhere  other  than  B1 
since  the  A  matrix  in  Bl  is  changed  to  an  identity  matrix  by  this  subroutine.  Instead 
of  testing  for  a  zero  determinant,  the  subroutine  uses  the  error  identification  subroutine 
9900  to  determine  if  the  A  matrix  is  not  invertablc.  Variables  K2-K9  are  used  to  avoid 
conflict  with  other  counters. 


8000  'STORE  Al(l,,) 

8010  R( 3  )=R( 1 1  :C<  3  )=C< 1 ) : FORIl=lTOR<  3 ): FORI2*lTOC< 3 i 
8012  Al( 3,11,12  )=A1( 1,11,12)  :NEXTI2:NEXTI1: RETURN 
8200  'RETRIEVE  THE  STORED  MATRIX 

8210  R( MI  )=R<  3  )  :C( MI  )=C( 3  ) : F0RI1=1T0R( MI  ) : F0RI2=1T0CI MI ) 
8212  A1(MI,I1, 12  )=A1< 3,11,12  )  :NEXTI2 :NEXTI1: RETURN 


9800  ' Simultaneous  Linear  Equation  Subroutine:  Ax=b 

9802  DIM  BKK9+1,K9*2),B2<K9J: 'B1  =  A  matrix}  B2  =  b  vector 

9805  'Input  from  SLEIN.DO;  Set  MAXFILES  in  main  program 

9806  OPEN"SLEIN" F0RINPUT AS9 

9807  F0RK8= 1T0K9 : F0RK7=1T0K9 : INPUT#9 ,B1< K8 ,K7 ) : NEXTK 7 : NEXTK8 

9808  F0RK8= 1TOK9: INPUT#9 ,B2  <  K8  ) : NEXTK8 
9815  'Invert  Matrix  A 

9820  F0RK7=K9+1T02*K9:F0RK8=1T0K9 
9822  IFK7=K8+K9THENB1< K8,K7  )=1ELSEB1( K8.K7  )=0 
9825  NEXTK8-.NEXTK7 
9830  F0RK7=1T0K9 

9835  IFBK K7,K7 )*SGN( BJ(K7>K7)  )<1E -8THENGOSUB9910 
9840  F0RK6=1T02*K9:BKK7,K6)=B1(K7,K6)/BKK7,K7):NEXTK6 
9842  IFK7=K9THEN9865 

9845  F0RK8=K7+1T0K9 :  IFBK  K8  ,K7  )=0THEN9860 

9855  F0RK6=K7T02*K9:B1( K8 ,K6 )=B1( K8 ,K6 J-t Bit K8  ,K7 )*B1< K7 ,K6 )) : NEXTK6 
9860  NEXTK8 : NEXTK  7 
9865  F0RK7=K9T02STEP-1 

9870  F0RK8=K7-1T01STEP-1 :  IFBK  K8.K7  }=0THEN9885 

9880  F0RK6=1T02*K9:BKK8,K6)=BKK8,K6)-<B1<K8,K7)*BKK7,K6)):NEXTK6 

9885  NEXTK8 : NEXTK  7 

9890  'Ffcilt  A  Inverse  by  b 

9892  PRINT'***  Sim  Lin  Eq  Solution:  X(i)  =" 

9894  F0RK7=1T0K9:B1(K7,1)=1:F0RK8=K9+1T02*K9:F0RK6=1T0K9 
9896  B1(K7>1  )=B1(K7>1  )+Bl(  K7>K8  )*B2(  K6  KNEXTK6  :NEXTK8 
9898  PRINTUSING"####. ###"  jBK K7 ,1 )  :NEXTK7:  PRINT"" : RETURN 
9900  'Error  Routine 

9903  IFERL>9700AN0ERR=11THENPRINT"!!  'ERROR:  Matrix  Is  Not  Inver  table  !!!*':  END 
9905  PRINT"Error  Code" } ERR} "In  Line" }ERL :EN0 
9910  'SWITCH  ROWS 

9915  F0RK5=K7+1T0K9:IFB1(K5,K7)*SGN(BKK5,K7)  X1E-8THEN9940 
9920  F0RK4=1T0K9*2:K3=BKK7,K4):BKK7,K4)=B1IK5,K4) 

9930  Bl(  K5  r  K4 ) =K3 : NEXTK4 : RETURN 

9940  NEXTK5 : PRINT"Error :  Matrix  Not  Invertable" : END 


VIII.  NUMERICAL  DOUBLE  INTEGRATION  PROGRAM 


A.  GENERAL 

This  program  numerically  integrates  functions  of  one  or  two  variables  using 
Simpson's  Rule  with  a  Romberg  extrapolation  to  improve  accuracy.  The  Romberg 
extrapolation  is  described  in  [Ref.  12:pp. 250-276].  Using  the  Romberg  extrapolations 
allows  the  operator  to  specify  an  acceptable  error.  The  program  conducts 
extrapolations  until  the  error  of  the  numerical  estimate  is  below  that  specified 
tolerance. 

The  operator  may  interactively  change  the  function  being  integrated,  the  limits  of 
integration,  or  the  Romberg  tolerance.  This  program  is  also  written  as  a  subroutine. 
However,  in  the  subroutine  the  function  being  integrated,  the  limits  of  integration,  and 
the  tolerance  may  not  be  interactively  changed. 

B.  INPUT 

All  input  is  entered  from  the  keyboard  of  the  M100.  When  the  program  begins, 
a  menu  appears  which  allows  the  operator  to  select  whether  the  function  to  be 
integrated,  the  limits  of  integration,  or  the  Romberg  tolerance  will  be  changed. 

When  the  operator  selects  an  input  to  be  changed,  the  program  calls  the  edit 
function  for  the  applicable  lines  in  the  program.  The  edit  function  terminates  the 
running  of  the  program.  The  operator  should  change  only  the  right  hand  side  of  the 
input  equations.  After  changing  the  lines  required,  the  operator  will  hit  the  F8  button 
on  the  M100.  This  puts  the  M100  back  in  the  BASIC  mode.  The  operator  must  then 
enter  RUN  or  hit  the  F4  button  to  run  the  program  again.  If  the  operator  wants  to 
change  another  input,  he  should  select  another  input  from  the  menu  and  repeat  the 
process. 

1.  Changing  The  Function,  f(x,y),  To  Be  Integrated 

Enter  zero  at  the  main  menu.  When  lines  1285-1288  appear,  change  the  right 
hand  side  of  the  equation  on  line  1286.  If  the  equation  is  too  long  for  one  line,  then: 

•  Calculate  a  partial  function  value  on  line  1286,  assigning  it  to  F. 

•  Add  a  line  1287  assigning  the  final  function  value  to  F  and  including  the  partial 

function  value  from  line  1286  in  the  right  hand  side  of  the  equation  on  line  1287. 
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For  example,  if  flXy)  =  (x  +y  +7)  *  eA  ^  ,  then  the  function  might  be  broken 
down  as  indicated  in  Figure  8.1. 


Figure  8.1  Example  Of  Function  To  Be  Integrated. 

Up  to  two  independent  variables,  X  and  Y,  may  be  used  in  the  equation.  The  operator 
must  ensure  that  f(x,y)  is  formulated  with  X  as  a  variable  for  which  constant  limits  of 
integration  can  be  specified.  After  f{x,y)  is  entered,  depress  F8,  then  F4  to  return  to 
the  main  menu. 

2.  Changing  The  Limits  Of  Integration 

Enter  one  at  the  main  menu.  When  lines  1291-1298  appear,  change  the  right 
hand  side  of  the  equations  on  lines  1293-1298  as  desired.  The  upper  and  lower  limits 
of  integration  for  X,  XUPPER  and  XLOWER,  must  be  constants.  The  upper  and 
lower  limits  of  integration  for  Y,  YUPPER  and  YLOWER,  may  be  constants  or  given 
in  terms  of  X.  Do  not  alter  the  return  statement  at  line  1295.  After  the  limits  of 
integration  are  entered,  depress  F8,  then  F4  to  return  to  the  main  menu. 

3.  Changing  the  Romberg  Tolerance 

The  operator  should  enter  two  from  the  main  menu  and  enter  the  new 
tolerance  when  prompted. 

4.  Using  The  Program  For  Single  Integration 

Although  the  program  is  written  for  double  integration,  single  integration  may 
be  calculated  using  the  following  steps. 

•  Set  the  function  to  be  integrated  equal  to  one,  i.e.  line  1286  will  be  F  =  1. 

•  In  lines  1296-1297  set  YUPPER-  fl[X)  and  YLOWER  =  0. 

•  In  lines  1293-1294  set  XUPPER  and  XLOWER  as  the  constant  limits  of  X 

between  which  the  function  YUPPER=fj[X)  is  to  be  integrated. 

2  7 

For  example,  for  J  dx,  the  corresponding  limits  of  integration  in  lines  1293-1297 
would  be  as  indicated  in  Figure  8.2. 
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1293  XUPPER=2 

1294  XLOWER=0 

1295  RETURN 

1296  YUPPER=Xa  2 

1297  YL0HER=0 


Figure  8.2  Example  of  Limits  Of  Integration. 

C.  OUTPUT 

The  estimated  value  of  the  integral  is  printed  to  the  screen  with  its  tolerance 
error.  If  the  program  generated  a  tolerance  error  that  was  less  than  the  tolerance 
specified  in  the  input,  then  that  tolerance  is  printed.  If  the  program  could  not  generate 
an  estimate  within  the  specified  tolerance,  then  a  message  to  that  effect  is  printed  to 
the  screen. 

D.  EXPLANATION  OF  VARIABLES 

•  A2(6,6)  is  the  matrix  holding  Romberg  extrapolation  values. 

•  DX  and  DY  are  the  widths  of  intervals  (XU-XL)/N  and  (YU-YL)/N  respectively. 

•  F  is  the  value  of  Rx,y)  to  be  integrated  at  a  particular  point. 

•  J1  through  J9  are  loop  counters. 

•  N  is  the  number  of  intervals  into  which  the  distances  XU-XL  and  YU-YL  are 
divided. 

•  SS  is  the  Simpson's  Rule  sum  specified  in  equation  8.1.  In  equation  8.1 
fj_  y|x=  x)  YL  —  ^  YU,  i=  1,2,3,. ..,n+  1.  n  is  the  number  of  intervals  into 
which  the  distance  YU-YL  has  been  divided,  n  must  be  a  positive,  even  integer. 

Simpsons's  Rule  Sum  =  fj +4f2  +  2fy  +  4f^4-  2fj +  ,...,  + 4fn  +  fn  +  j  (eqn  8.1) 

•  TL  is  the  user  specified  tolerance. 

•  XLOWER  or  XL  is  the  lower  limit  of  integration  for  X. 

•  XUPPER  or  XU  is  the  upper  limit  of  integration  for  X. 

•  YLOWER  or  YL  is  the  lower  limit  of  integration  for  Y. 

•  YUPPER  or  YU  is  the  upper  limit  of  integration  for  Y. 

•  Z9  is  a  selection  variable. 
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E.  EXPLANATION  OF  PROGRAM  COMPONENTS 

A  complete  program  listing  is  located  at  Appendix  F. 

1.  Initialization,  Figure  8.3 


Figure  8.3  Initialization  Section. 

Line  1201  dimensions  the  matrix  holding  the  Romberg  extrapolations  and  sets 
the  default  tolerance  to  .001. 

2.  Option  Selection,  Figure  8.4 

1205  CLS: PRINT"": PRINT"  *«  Double  Integration  **" 

1206  PRINT"  Romberg  Algorithm" 

1210  PRINT"0=Edit  Function  To  Be  Integrated." 

1211  PRINT"l=Edit  Limits  Of  Integration." 

1213  PRINT"2=Edit  Tolerancej  Current  Tol.="jTL 

1215  RRINT"3=Calculate  Integral  " :INPUT"Enter  0,  1,  2,  or  3:"»Z9 

1216  IFZ9=0THENE0IT1285-1288 

1217  IFZ9=1THENEDIT1291-1298 

1218  IFZ9=2THENPRINT"":INPUT"Tolerance="»TL:G0T01205 

Figure  8.4  Option  Selection  Section. 

Lines  1205-1218  print  a  menu  which  permits  the  operator  to  change  f(x,y),  the 
limits  of  integration,  or  the  tolerance,  or  to  calculate  the  integral.  If  f(x,y)  or  the  limits 
of  integration  are  to  be  changed,  then  lines  1216  or  1217  activate  the  editor  for  the 
appropriate  program  lines.  The  editor  terminates  program  execution,  thereby  requiring 
that  the  program  be  executed  after  editing.  If  the  tolerance  is  to  be  changed,  then  line 
1218  prompts  the  operator  to  update  TL  and  redisplays  the  menu. 

3.  Integration  Calculation,  Figure  8.5 

Line  1220  clears  the  screen  during  the  calculation  and  prints  an  admonition  to 
be  patient  while  the  calculation  occures.  Line  1230  sets  the  initial  number  of  intervals 
to  two,  calls  subroutine  1293,  which  calculates  the  interval  width,  DX.  Line  1240 


1220  CLS: PRINT"": PRINT"  !'B«  Patient! !":PRINT"" 

1230  N=2:GOSUB1293:DX=<XU-XU/2 

1240  F0RJ9=1T06 : DX=DX/2 : N=N*2 

1242  X=XU:GOSUB1296:GOSUB1280:A2« J9,1)=SS*DY 

1245  X=XL : 60SUB1296 : GOSUB1280 : AZ l J9 , 1 > = A2 ( J9 , 1 ) +SS*DY 

1250  FORJ8=2TON:X=X+DX:GOSUB1296 : GOSUB1Z80 

1251  A2( J9,1)=AZ( J9,1I+2*SS*0Y:NEXTJ8 

1252  A2(J9,1)=A2( J9,l)#0X/3 


Figure  8.5  Integration  Calculation. 


starts  a  loop  in  which  the  Simpson's  Rule  intervals  are  halved  at  each  iteration.  That 
is,  in  the  first  iteration  YU-YL  and  XU-XL  are  divided  into  four  intervals,  in  the 
second  iteration  y  are  divided  into  eight  intervals,  and  so  on  for  six  iterations.  Line 
!  1240  cuts  the  interval  for  X,  DX,  in  half  and  doubles  the  number  of  intervals,  N. 

i 

Lines  1242-1245  call  the  subroutines  which  compute  the  Simpson  Rule  sums, 
SS,  at  the  upper  and  lower  bounds  of  X.  These  sums  are  multiplied  by  their  respective 
interval  widths,  DY,  and  added  together  into  A2(J9,1).  Lines  1250-1251  calculate  the 
!  same  summation  for  values  of  X  between  XL  and  XU  at  intervals  DX  and  and  add  the 

5  sums  to  A2(J9,1).  Line  1252  multiplies  A2(J9,1)  by  DX/3  to  complete  the  Simpson's 

Rule  approximation  of  f(x,y)  dydx. 

4.  Romberg  Extrapolation,  Figure  8.6 


1255  IFJ9=1THENNEXTJ9 

1260  F0RJ8=1T0J9-1 

1261  A2(  J9,J8+1  )=A2(  J9, J8  )♦<  ( A2(  J9, J8  )-A2(  J9-1,  J8  J )/( 4A J8-1 ) ) -.NEXTJ8 


Figure  8.6  Romberg  Extrapolation. 


Because  the  Romberg  extrapolations  require  two  numerical  approximations, 
line  1255  skips  the  extrapolation  section  after  the  first  iteration,  i.e.  when  J9=  1.  Lines 
1260-1261  conduct  the  Romberg  extrapolation  as  described  in  the  section  on  Romberg 
extrapolation  in  [Ref.  12:pp.  250-276]. 


5.  Termination  and  Output,  Figure  8.7 
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T1=A2( J9»J9  )-A2( J9, J9-1 ) : IFSGNI T1 )»T1-TL>0THENNEXTJ9ELSE1264 
PRINT"Toloranca  of"jTL>"not  mot  after  five  extrapolations” 

IN=A2< J9,J9 ) 

PRINT  "Integral  ="  J :  PRINTUSING"###### .  #####••  jIN 

PRINT"  Actual  Toleranco="» : PRINTUSING"##. ####">T1*SGN( Tl ) 

SOUND 1567,10: SOUND 1244,10: SOUND 1046 , 10 : SOUND  783,20 
SOUND 1046 , 10 : S0UND783 ,40 
INPUT"Hit  Enter  To  Continue:" jZ9:GOT01205 
F0RJ7=1T06 : F0RJ6=1T0J7 : PRINTUSING"## . ###" j A2( J7, J6 ) > 

NEXT J6 : PRINT NEXT J  7 : INPUTZ9 : RETURN 


Figure  8.7  Program  Termination  and  Output. 


Line  1262  finds  the  difference,  Tl,  between  the  last  two  Romberg 
extrapolations  and  compares  that  difference  to  the  user  specified  tolerance,  TL.  If  the 
difference  is  greater  than  the  tolerance,  then  the  extrapolation  is  not  accurate  enough, 
another  numerical  integration  is  conducted  with  DX  and  DY  halved,  and  another 
extrapolation  is  made.  Otherwise,  the  integral  is  within  tolerance  and  the  program 
branches  to  line  1264  to  begin  the  output  sequence.  If  the  integral  is  not  within 
tolerance  after  six  iterations,  iterations,16  the  program  terminates.  Line  1263  prints  a 
message  to  the  screen  indicating  that  the  tolerance  has  not  been  met.  Line  1264 
assigns  the  final  value  of  the  integral  to  IN.  Lines  1265  and  1266  print  the  value  of  the 
integral  and  the  actual  tolerance,17  to  the  screen.  Lines  1267  and  1268  play  a  short 
tune  to  cue  the  operator  that  the  calculation  has  finished.  Line  1269  holds  the  results 
on  the  screen  until  the  operator  hits  ENTER,  cycling  the  program  back  to  the  main 
menu  at  line  1205. 

6.  Diagnostic  Subroutine,  Figure  8.8 

Lines  1275-1276  print  the  A2  matrix.  This  subroutine  can  be  called  in  the 
middle  of  a  calculation  to  check  how  far  the  calculation  has  progressed.  To  call  the 
subroutine  in  the  middle  of  a  calculation: 

•  Hit  SHIFT  and  BREAK  together  to  stop  the  calculation. 

•  Enter  GOSUB 1275. 


16That  is,  after  distances  XU-XL  and  YU-YL  have  been  broken  into  128 
intervals. 


17Actual  toleranc  may  be  less  than  the  user  specified  tolerance. 


1275  F0RJ7=1T06 : F0RJ6=1T0J7 : PRINTUSING"## . ###" J  A2  <  J7 , J6 ) 

1276  NEXT J6 : PRINT"" : NEXT J7 : INPUTZ9: RETURN 


Figure  8.8  Diagnostic  Subroutine,  Prints  Matrix  A2.. 

•  After  viewing  the  matrix,  hit  ENTER  to  continue  the  program. 

7.  Simpson's  Rule  Summation,  Figure  8.9 

1280  REM  Simpson’s  Rule  Sum 

1281  Y=YU:G0SUB1285:SS=F:Y=YL:G0SUB1285:SS=SS+F 

1282  F0RJ5=2T0( N/2  ) : Y=Y+DY : G0SUB1 285 : SS=SS+4*F : Y  =Y+DY : G0SUB1285 

1283  SS=SS+2*F : NEXTJ5 : Y=Y  +0Y : G0SUB1285 : SS=SS+4*F : RETURN 


Figure  8.9  Simpson's  Rule  Summation  Subroutine. 

Lines  1281-1283  calculate  SS=  £  f j  +  4^  +  2 fj  +  4f^  +  2f^  +  ,...,  4fn  + 
fn+l  where  ^  =  f(xy|x=X)  YL£yj5YU,  i=  1,2,3,. ..,n+ 1.  n  is  the  number  of 
intervals  into  which  the  distance  YU-YL  has  been  divided. 

8.  F(x,y)  to  be  integrated,  Figure  8.10 


1285  ’  f(x,yl  to  be  integrated: 

1286  F=1 

1288  ’X  i  Y=independent  variables.  Hit  F8,  then  F4  When  Done. 

1289  RETURN 

Figure  8.10  Subroutine  To  Calculate  f(x,y). 

The  function  to  be  integrated,  H[x,y)  is  at  line  1286. 18  Lines  1285  and  1288  are 
comments  printed  to  the  screen  during  editing  to  assist  the  operator. 

18An  additional  line,  1287,  may  be  added  if  the  function  is  too  long  for  one  line. 


9.  Limits  Of  Integration,  Figure  8.11 


1290  'Limits  of  Integration: 

1291  1 XLOHER/XUPPER  are  constants. 

1292  'YUPPER  S  YLOWER  may  be  constants  or  given  in  terms  of  X. 

1293  XUPPER=1. 5707963 

1294  XL0WER=0 

1295  RETURN 

1296  YUPPER=SIN<  X ) 

1297  YLOWER=0 

1298  ‘Hit  F8>  Then  F4  When  Done 

1299  DY=(YU-YL)/<N+1):  RETURN 


Figure  8.1 1  Limits  Of  Integration  Subroutines. 

Upper  and  lower  limits  of  integration  for  X  are  entered  at  lines  1293  and  1294 
respectively  and  must  be  constants.  Limits  of  integration  for  Y  are  entered  at  lines 
1296-1297  and  may  be  either  constants  or  functions  of  X.  Line  1299  updates  DY. 
Lines  1290-1292  and  1298  are  comments  to  assist  the  operator  during  editing. 

F.  INTEGRATION  SUBROUTINE 

The  numerical  integration  program  described  above  is  adapted  in  Figure  8.12  for 
use  as  a  subroutine.  In  the  subroutine  neither  f(x,y),  the  limits  of  integration,  nor  the 
tolerance  can  not  be  edited  during  program  execution.  All  comment  lines  to  facilitate 
editing  have  been  removed.  The  subroutine  returns  IN  as  the  numerical  approximation 
of  the  integral  but  does  not  print  IN.  The  operator  must  dimension  A2(6,6)  with  the 
other  arrays  in  the  main  program  and  delete  line  1201  in  the  subroutine. 


1200  'Numberical  Integration  Subroutine: Steven  H.  Cary: 24  Apr  86 

1201  DIMA2I 6>6 ) 

1202  TL=.001 

1220  CLS: PRINT"": PRINT"  ! 'Calculating  An  Integral! !":PRINT"" 

1250  N=2:GOSUB1293:OX=<XU-XL)/2 

1240  F0RJ9=1T06 : 0X=DX/2 : N=N*2 

1242  X=XU:GOSUB1296:GOSUB1280:A2f J9,l  )=SS*DY 

1245  X=XL : G0SUB1296 : G0SUB1280 : A2 ( J9 , 1 >=A2 ( J9 , 1 )  +SS*DY 

1250  FOR J8= 2T0N : X=X+0X : GOSUB 1296: GOSUB 1 280 

1251  A2( J9,1)=A2( J9,l  )+2*SS*DY :NEXTJ8 

1252  A2< J9>1  )=A2( J9»l )*DX/3 
1255  IFJ9=1THENNEXTJ9 

1260  F0RJ8=1T0J9-1 

1262  A2C  J9 , J8+1  )=A2(  J9 ,  J8  )♦(  ( A2(  J9 ,  J8  >-A2(  J9-1 ,  J8  ))/< 4A  J8-1 ) )  :NEXTJ8 

1263  T1=A2<  J9,J9)-A2<  J9.J9-1  ):IFSGN(  T1  )»T1-TL>0THEWEXTJ9ELSE1266 

1264  PRINT"Tolerance  of" iTL ("not  met  after  five  extrapolations" 

1266  IN=A2(  J9.J9): RETURN 

1275  F0RJ7=1T06 : F0RJ6=1T0J7: PRINTUSING"## . ###" J A2<  J7,J6  ) i 

1276  NEXTJ6 : PRINT"" :NEXTJ7: INPUTZ9 : RETURN 

1280  'Simpson’s  Rule  Sum 

1281  Y =YU : GOSUB 1 285 :SS=F : Y=YL :G0SUB1285:SS=SS+F 

1282  F0RJ5=2T0( N/2  ) : Y=Y+DY : G0SUB1285 : SS=SS+4*F : Y=Y+DY : GOSUB1285 

1283  SS=SS+2*F : NEXT J5 : Y=Y+DY :G0SUB1285:SS=SS+4*F: RETURN 
1286  F=1 

1289  RETURN 

1293  XUPPER=1 

1294  XLOWER=0 

1295  RETURN 

1296  YUPPER=X 

1297  YLOWER=0 

1299  DY=<YU-YL)/<N+1):  RETURN 


APPENDIX  A 

DETECTION  SIMULATION  PROGRAM  LISTING 


A  complete  listing  of  the  Detection  Simulation  Program  is  as  follows. 


100  CLS :PRINT"": PRINT"  DETECTION  SIMULATION" : FORI=1T0400:NEXTI 

110  'Input/Initialization 
115  0PEN”DSIN"F0RINPUTAS1 

120  INPirntl ,NS,NP ,S1,S2 ,RH ,F1 : V1=S1*S1: V2=S2*S2 

121  DIMX(NS,5+NP),A2<6,6),T113) 

125  FORI 1 = 1T05+NP : INPUT#1 »X( 1,11): NEXTI1 

126  IFNS=1THEN140 

127  IFF1=1THEN132 

128  FORI l=2TONS: FORI 2=1T05+NP:INPUT#1>X( II >12  )  :NEXTI2:NEXTI1 :GOT0140 
132  FORI1=2TONS:FORI2=1T05:INPUT#1,X(I1,I2):NEXTI2:IFNP=OTHEN135 

134  FORI 2=6T05+NP :  XI 1 1 , 12 ) =X< 1 ,1 2  ) : NEXTI 2 

135  NEXTI1 

140  RF=SQR(1-RHa2) 

150  DF$="EXP(  -( ( XT-XS  lA2+<  YT-YS  )A2  )/(  2«XU2»6  )A2  ) )" 

200  '***  Simulation  Selection  Section  *** 

201  CLS:PRINT"Is  the  Detection  function:" 

203  PRINT"  1.  Deterministic" : PRINT"  2.  Probabilistic" 

205  INPUT"Enter  1  or  2:"jFl 

210  CLS:PRINT"Are  Sensor  Locations:" 

212  PRINT"  1.  Always  At  Aim  Point" : PRINT"  2.  Distributed  BVN  Around  Aim  Point" 

214  INPUT"Enter  1  or  2:"}F2 

215  IF( F1=20RF2=2 )THENF3=1 :GOTO230 

220  CLS: PRINT"" :PRINT"Is  the  Calculation:" 

222  PRINT"  l=Monte  Carlo  Simulation" : PRINT"  2=Numerical  Approximation” 

224  INPUT"Enter  1  or  2">F3 

230  TIME$="00:00:00":IF  F1=1THENGOSUB300ELSEGOSUB500 
250  GOTO200 

300  'Deterministic  Sensor  Subroutine 

305  I FF3= 1THENGOSUB310E  LSEGOSUB350 

306  RETURN 

310  'Monte  Carlo  of  Deterministic  Sensor 
315  GOSU8900 

320  PD=0 : F0RJ1=1T0NR : PRINT . 241 , "Repe t i t ion: " > J1 : GOSUB600 : F0RJ2=1T0NS 

323  IFF2=2THENXS=X( J2 ,1 ) : YS=X( J2 ,2 ) :G0T0325 

324  GOSUB612 

325  T1=SQR(  (XS-XT  )A2+<  YS-YT  lA2  ) 

330  IFT1<=X( J2 ,6  )THENPD=PD+1 :G0T0335 

332  IFT1>=X( J2 , 7  )ANDT1<=X( J2 ,8  )THENPD=PD+1 : G0T0335 

334  NEXTJ2 

335  NEXTJ1 : PD=PD/NR : GOSUB950 : RETURN 
350  'Numeric/Deterministic  Subroutine 

355  P0=0 : F0RJ2=1T0NS: H=X( J2 ,6 ) :GOSU81200 : PD=P0+IN 

356  H=X( J2,8):GOSUB1200:PD=P0+IN 

357  H=X( J2>7 ) : GOSUB1200 : PD=P0-IN:NEXTJ2 
360  GOSUB950: RETURN 

500  'Probabilistic  Detection  Function 

502  CLS:PRINT"Qefault  Detection  Function  Is  Carleton." 

503  GOSUB1300:GOSUB900 

520  P0=0 : F0RJ1=1T0NR : PRINT . 241 , "Repetition : " > J1 : GOSUB600 : F0RJ2=1T0NS 

521  IFF2=2THENXS=X( J2 ,1 1 : YS=X( J2 ,2  I  :G0T0523 

522  GOSUB612 

523  GOSUB1410 : IFRND( 1  )<=DFTHENPD=PD  +  1 :G0T0526 

524  NEXTJ2 

526  NEXTJ1:P0=PD/NR 
530  GOSUB950: RETURN 
600  '***Generate  BVN  RV»*» 

602  U1~RN0( 1 )  :U2=RN0( 1  ):TE=SQR( -2*L0G(U1  )  ) 
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604  XT=TE*C0S(  6 . 2831853*112  ) :  YT=RH*XT+RF*TE*SIN<  6 . 2831853*112  ) 

606  XT=XT*S1:YT=YT*S2: RETURN 

612  U1=RND( 1 )  ;U2=RND< 1 )  :TE=SQR( -2*L0G( U1 ) ) 

614  XS=TE*C0S( 6 . 2831853*U2 ) : YS=X(  J2 ,5 )*XT+< 1-Xt  J2 ,5  )*2 ) A . 5*TE*SIN( 6 . 2831853*U2 ) 
616  XS=X( J2 , 1  )+XS*X< J2 ,3 ) : YS=X<  J2 , 2  )+YS*X<  J2 ,4 ) : RETURN 

900  CLS:INPUT"Enter  number  of  repetitions  for  Monte  Carlo  Simulation: " >NR 
905  RETURN 

910  INPUT"Hit  ENTER  to  Continue” )Z1 : RETURN 

950  ‘Print  output 

951  SOUND 156 7 , 10 : SOUND 1244 , 10 : SOUND 1046  > 10 : SOUND  783,20 

952  SOUND 1046 ,10 : SOUND 783 ,40 

953  CLS: PRINT"" :PRINT"Calculation  Time  (HH/IWSS)  =  “ >TIME$ :IFF3=2THEN960 

954  PRINT"Select  Alpha  for  Confidence  interval:" 

955  INPUT"  Choices  =  .1,  .05,  .01:")AL 

956  IFAL= . 1THENAL=1 . 645 : GOT0960 

957  IFAL=.05THENAL=1. 96:GOT0960 

958  IFAL= . 01THENAL=2 . 575 : GOT0960 

959  G0T0954 

960  PRINT"***  Estimate  of  P( Detection)  =  “ > : PRINTUSING"  #.#####" )PD 

961  IFF3=1THEN965 

962  PRINT"No  Confidence  Interval  For  Numerical  Approximations" 

963  GOT0970 

965  PRINT "Confidence  Interval:  ("j 

966  TE=AL*SQRl PD*( 1-PD  )/NR )  :LL=PD-TE :UL=PD+TE :IFUL>1THENUL=1 

967  IFLL<0THENLL=0 

968  PRINT  USING1'))# .  ###))))"  jLL>UL) :  PRINT"  )":  G0SUB910 
970  ‘Confetti  Approximation 

972  PRINT"" :INPUT“Confetti  approximation?  0=No,  l=Yes:“»Z9:IFZ9=0THENRETURN 
974  CLS:INPUT”Enter  TOTAL  lethal  area  for  ALL  sensors  in  the  pattem:“»NA 

976  TE=NA/( 6 . 283185*S1*S2  )  :TE  =  l-( l+SQRI 2*TE  )  )*EXP( -SQR< 2*TE  ) ) 

977  PRINT“**Confetti  Approximation  =  “ »TE :GOSUB910: RETURN 

1200  'Numberical  Integration  Subroutine 

1201  D1=6.2831853*S1*S2*RF 

1202  TL=.001 

1220  CLS: PRINT"": PRINT"  ! 'Calculating  An  Integral! !": PRINT"" 

1230  N=2 :GOSUB1293 : DY=l YU-YL )/2 

1240  F0RJ9=1T06:DY=DY/2:N=N*2 

1242  Y=YU:GOSUB1296:GOSUB1280:A2(J9,1)=TS*DX 

1245  Y=YL : G0SUB1296 : GOSUB1280 : A2( J9 , 1  )=A2( J9 ,1  HTS*DX 

1250  F0RJ8=2T0N : Y  =Y*0Y : GOSUB1296 : GOSUB1280 

1251  A2( J9,l  )=A2( J9,l  )+2*TS*0X : NEXTJ8 

1252  A2<  J9,l  )=A2(  J9,l  )*0Y/2 
1255  IFJ9=1THENNEXTJ9 

1260  FORJ8=lTOJ9-l 

1262  A2(  J9,J8+I  )=A2(  J9,J8  )♦(  ( A2(  J9,  J8  )-A2(  J9-1  ,J8  )  )/(4AJ8-l )  ):NEXTJ8 

1263  T1=A2( J9,J9 )-A2( J9 ,  J9-1 ) : IFSGNt T1  )*T1-TL>0THENNEXTJ9ELSE 1266 

1264  PRINT"Tolerance  of">TL)"not  met  after  five  extrapolations" 

1266  IN=A2(J9,J9):  RETURN 

1275  F0RJ7=1T06 :  FORJ6=lTOJ7 :  PRINTUSING")))) .  ###"  )  A2 (  J7 ,  J6  ) ) 

1276  NEXTJ6 : PR  lNT“" :NEXTJ7 : INPUTZ9 : RETURN 

1280  REM  Trapezoidal  Rule  Sum 

1281  X=XU:GOSUB1286:TS=F:X=XL:GOSUB1286:TS=TS*F 

1282  F0RJ5=2T0N-1:X=X+DX:G0SUB1Z86 : TS=TS*F : NEXT J5: RETURN 

1285  'F(X,Y)  to  be  integrated: 

1286  F=XA2/V1-2*RH*X*Y/S1/S2*YA  2/V2 

1287  F=(EXPt-F/2/RFA  2  n/01:  RETURN 
1290  'Limits  of  Integration: 

1293  YU=XI  J2,2  HH: YL=X(  J2 ,2  )-H:RETURN 

1296  T3=SQR(Ha2-(Y-X(  J2,2  (  lA  2  )  :XU=X<  J2 ,1  HT3:XL=X(  JZ  ,1  )-T3:DX=(XU-XL)/N 

1297  RETURN 

1300  PRINT"  -Detection  Fn  (OF)  in  terms  of  XT,  YT," 

1302  PRINT"  and  Parameters  XS,  YS,  and  X(J2,6) _ X(J2,5+NP):“ 

1304  PRINT"  **  OF  =  "»DF* 

1306  PRINT"Hit  ENTER  For  No  Change  or  Enter  New. . . “ :INPUT"  DF  =  ">DF$ 

1307  RETURN 

1400  'Tokenize  OF 

1410  B$="DF=“+DFS*CHR»(0) 

1450  'Tokenize/cxecute  89 

1451  B0=VARPTR(B$ )  :B1=PEEK< B0*1 H256*PEEKl  B0  +  2 ) :CALL1606 ,0,B1 
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CALL 2499, Oi 631 05 .-RETURN 


APPENDIX  B 

KALMAN  FILTER  PROGRAM  LISTING 


A  complete  listing  of  the  Kalman  Filter  Program  is  as  follows. 


100  CLS:PRINT"*****KALMAN  FILTER*****" : PRINT"  Input  Data  Baing  Read" 

110  0PEN"KALIN"F0RINPUTAS1 :ONERRORGOTO9900 
120  INPUT81 , NX , NZ : I FNX<NZTHENMO =NZE  LSEMD =NX 

125  DIMPH(NX,NX),MH(NX),Q<NX,NX),H<NZ,NX),MV(NZ),R1<NZ,NZ),MU(NX),SG<NX,NX) 

126  OIMCK  MD  *MD  )  ,C2(  MD  ,MD  )  »K(  NX  >NZ ) 

127  DIM81(NZ+1,NZ*2) 

130  F0RI1=1T0NX:  F0RI2=lT0NX:INPLfT#l,PH(  II  >12  )  :NEXTI2 :NEXTI1 
132  FORI  1=1T0NX:  INPUT It  1 ,MH( II )  :NEXTI1 

134  FORIl=lTONX:FORI2=lTONX:INPUTttl,Q< 11,12  )  :NEXTI2:NEXTI1 
136  F0RI1=1T0NZ:  F0RI2=1T0NX:  INPUT  tfl ,H(  II »I2 ) :NEXTI2 :NEXTI1 
138  FORIl=lTONZ:INPUTttl,MVIIl):NEXTIl 

140  F0RI1=1T0NZ:  F0RI2=1T0NZ:  INPUTItl,Rl(  II  >12  )  :NEXTI2:NEXTI1 
142  F0RI1=1T0NX: INPUT81 ,MUI II )  :NEXTI1 

144  F0RI1=1T0NX : FORI  2=1T0NX :  INPUTIfl,  SGI  11,12  )  :NEXTI2 :NEXTI1 

145  CC=0:CLS: PRINT “Initial  SG  As  Input  Check:":GOSUB532 

150  CLS: PRINT"  *****MEASUREMENT  BLOCK*****" 

160  PRINT"Current  H  ":GOSUB540 

162  INPUT-Enter  New  H  ?  l=Yes,  0=No:")Z9:IFZ9=0THEN170 
165  'Enter  A  New  H 

167  F0RI1=1T0NZ:F0RI2=1T0NX 

168  PRINT"Enter  Row")Ilj",  Column" )I2 l"0f  H 

169  INPUTHI II ,12  )  :NEXTI2 : PRINT"" :NEXTI1 

170  'CALC  KALMAN  GAIN 

171  'MULT  SG  H  t,  INTO  Cl 

172  F0RI1= 1TONX : FORI 2 = 1TONZ : C 1 1 1 1 , 1 2  ) =0 : FORI 3 =1T0NX 

174  C1(I1,I2 J=< SGI  11,13 )*H(  12,13 ))+Cl< 11,12 J:NEXTI3:NEXTI2:NEXTI1 
180  'MULT  H  SG  H  t  ,  INTO  C2 

182  F0RI1=1T0NZ : F0RI2=1T0NZ:C21 II ,12  )=0 : F0RI3=1T0NX 

184  C2(  11,12  1  =  1  HI  11,13  I*C1(  13,12)  )+CZI  11,12  ):NEXTI3:NEXTI2:NEXTI1 

200  'ADD  R  INTO  C2 

202  F0RI1=1T0NZ: F0RI2=1T0NZ:C2< 11,12 )=C21 11,12 )+Rl(  11,12 ) 

203  NEXTI2 :NEXTI1 
210  'INVERT  C2 
215  GOSUB9800 

220  'MULT  Cl  C2  INTO  K 

222  F0RI1=1T0NX:F0RI2=1T0NZ:K( 11,12 )=0 : F0RI3=1T0NZ 

224  K(  11,12  )=(C1(  11,13  )*C2(  13,12  )  )+Kl  11,12  ) :NEXTI3 :NEXTI2:NEXTI1 

250  '*****UPDATE  MU-  TO  MU+**** 

251  'MULT  H  MU-  INTO  Cl 

252  F0RI1=1T0NZ:C1<  11,1 )=0 : F0RI3=1T0NX 

254  Cll  11,1  )=(  HI  11,13  )*MU(  13  )  1+Cll II ,1 ) :NEXTI3:NEXTI1 
260  'ADD  MV  ♦  H  MU- 

262  F0RI1=1T0NZ: Cll 11,1  )=C1( II ,1 )+MVI II ) :NEXTI1 
270  'INPUT  A  NEW  MEASUREMENT 

272  CC=CC+1 :CLS:PRINT”Measurement  tt")CC>":" 

273  FORIl=lTONZ:PRINT"Enter  Element" >11 »"Of  Measurement:") 

274  INPUTZI II ) :NEXTI1 

280  'SUBTRACT  Cl  FROM  Z,  INTO  Cl 

282  FORI 1=1T0NZ: Cll 11,1  l=ZI II  )-Cl< 11,1  ):NEXTI1 

290  'MULT  K  Cl  INTO  C2 

292  F0RI1=1T0NX:C2( 11,1  )=0 : F0RI3=1T0NZ 

294  C21 11,1  )  =  (K(  11,13  I*C1I  I3,1))+C2(I1,1 )  :NEXTI3  :NEXTI1 

300  'ADO  C2  ♦  MU-  TO  UPDATE  TO  MU+ 

302  F0RI1=1T0NX:MU( II )=C2( 11,1 )+MUI II  ):NEXTI1 

320  'MULT  K  H  *  SUBTR  FROM  I  ,  PUT  IN  Cl 

322  F0RI1=1T0NX:F0RI2=1T0NX:C1( 11,12 )=0:F0RI3=1T0NZ 

324  Cll  11,12  )=(KI  11,13  )*Hl  13,12  ))+Cl(  11,12  ):NEXTI3:C11 11,12  )=-Cl(  11,12  ) 
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326  NEXT I 2 : NEXTI 1 

328  F0RI1=1T0NX:C11 II >11 )=1+C1(  II ,11 ) : NEXTI 1 

350  ’MULT  LAST  RESULT  BY  SG  ,  INTO  C2 

352  FORI1=1TONX:FORI2=1TONX:C2II1,I2)=0:FORI3=1TONX 

354  C2<  11,12  )  =  (C1(  11,13  )*SG<  13,12  ) )+C2<  11,12  )  :NEXTI3  :NEXTI2:NEXTI1 

360  ’PUT  C2  INTO  SG 

362  F0RI1=1T0NX:F0RI2=1T0NX: SGI  11,12 )=C2( 11,12  ). -NEXTI 2: NEXTI 1 
375  CLS : PRINT"Kalman  Gain,  Kli.j)  After" 

377  PRINT "Measurement  8" >CC :GOSUB510 

380  CLS: PRINT"Estimate  Of  System  State,  MUI  i  )♦  After" 

382  PRINT "Measurement  8" >CC : GOSUB520 

385  CLS: PRINT"Estimate  Of  Covar ,  SGli,j)+  After" 

387  PRINT "Measurement  8" jCC:GOSUB530 

400  CLS:PRINT"*********MOV£MENT  BLOCK*******" 

410  ’Update  MUICO  +  to  MUICC  +  1)- 

420  ’MULT  PH  MU  ,  PUT  IN  Cl 

422  F0RI1=1T0NX:C1(I1,1  )=0 : F0RI3=1T0NX 

424  Cl<  11,1  )=t  PHI  11,13  )*MU<  13  )  I+C1I II ,1 ) :N£XTI3 :N£XTI1 

430  ‘ADD  Cl*  MW  ,  INTO  MU 

432  FORI 1=1T0NX: MUI II  )=C1(  II  ,1 )+MW(  II ): NEXTI 1  . 

440  ’**UP0ATE  SG  ** 

450  ’MULT  T  SG  ,  INTO  Cl 

452  F0RI1=1T0NX:F0RI2=1T0NX:C1< 11,12  1=0 : FORI3=1TONX 

454  Cll  11,12  )  =  <  PHI  11,13  )*SG<  13,12  )  )+Cl<  11,12  ):NEXTI3:NEXTI2:NEXTI1 

460  ’MULT  Cl  PH  t,  INTO  C2 

462  FORI1=1TONX:FORI2=1TONX:C2( 11,12  1=0 : F0RI3=1T0NX 

464  C2(I1,I2)  =  (C1(11,I3  )*PH(  12,13)  >♦02(11,12  ):  NEXTI  3 :  NEXTI  2 :  NEXTI  1 

470  ’ADD  C2  +  Q  =  SG 

472  FORI 1=1T0NX : FORI 2=1T0NX : SGI  II ,12  )=C2<  II ,12  )+qi  II ,12  )  :NEXTI2 :NEXTI1 
480  PRINT-Estimate  Of  System  State,  MUIII-" 

482  PRINT"Before  Measurement  8" |CC+1:GOSUB520 
485  CLS : PRINT"Estima te  Of  Covar,  SGII.J)-  Before" 

487  PRINT "Measurement  8">CC+1:GOSUB530 
490  GOT0160 

500  ’PRINTING  SUBROUTINES 
510  ’PRINT  KALMAN  GAIN,  K 

512  FORI1=1TONX : F0RI2=1T0NZ : PRINTUSING"888#88.88" )K(Il,I2)l : NEXTI 2 
514  PRINT"" : NEXTI 1 : INPUT “Hit  ENTER  To  Continue :" jZ9: RETURN 
520  'PRINT  MU 

522  F0RI1=1T0NX: PRINTUSING"888888 . 88" >MUI II ) \ : NEXTI 1 : PRINT"" 

524  INPUT"Hit  ENTER  To  Cont inue jZ9: RETURN 
530  ’PRINT  COVAR  MATRIX,  SG 

532  F0RI1=1T0NX: F0RI2=1T0NX : PRINTUSING"88888.88" >SG(I1,I2)> : NEXTI 2 
534  PRINT”" : NEXTI 1 : INPUT"Hi t  ENTER  To  Continue: " jZ9: RETURN 
540  ’PRINT  H 

542  F0RI1=1T0NZ: F0RI2=1T0NX : PRINTUSING"888888. 88" vHI 11,12)1 :NEXTI2 
544  PRINT’’"  : NEXTI  1 :  RETURN 
550  PRINT"  C2  MATRIX:" 

552  F0RI1=1T0A : F0RI2=1T0B : PRINT USING"888888 . 88" >C2I 11,12 )> :NEXTI2 
554  PRINT’"’ :  NEXTI  1 :  INPUT  "Hit  ENTER  To  Continue :” >Z9: RETURN 
9800  ’INVERT  C2 

9815  F0RI1=1T0NZ: F0RI2=1T0NZ:B1( 11,12 )=C2( II ,12  )  :NEXTI2 : NEXTI1 

9820  F0RI1=NZ+1T02*NZ:F0RI2=1T0NZ 

9822  IFI1=I2+NZTHENB1<  12,11 )=1ELSEB1(  12,11 )=0 

9825  NEXTI 2: NEXTI 1 

9830  F0RI1=1T0NZ 

9840  ML=1/B1<  11,11):  F0RI3=1T02*NZ:B11 11,13  )=Bll  11,13  )*ML :NEXTI3 
9842  IFI1=NZTHEN9865 

9845  F0RI2=I1+1T0NZ:IFB1< 12,11  )=0THEN9860 
9850  ML=-B1( 12,11) 

9855  F0RI3=I1T02*NZ:B1( 12,13 )=B1( 12,13 )  +  <ML*Bl( 11,13  )  ):NEXTI3 
9860  NEXTI2: NEXTI 1 
9865  F0RI1=NZT02STEP-1 

9870  F0RI2=I1-1T01STEP-1:IFB1( 12,11  )=0THEN9885 
9875  ML=-B1(  12,11) 

9880  F0RI3=1T02*NZ:B1( 12,13  )=B1( 12,13  )♦<  ML*B1(  11 ,13  )  ):NEXTI3 

9885  NEXTI 2: NEXTI 1 

9890  F0RI1=1T0NZ: F0RI2=1T0NZ 

9895  C21 12,11  )=B11 12,11+NZ  )  :NEXTI2:NEXTI1 
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9897  MI =1: RETURN 

9900  IFERL>9700AN0ERR=11THENPRINT"! ! 'ERROR:  C2  Is  Not  Invariable* ! ! " : END 
9905  PRINT"Error  Code" >ERR»"In  t  .ia"  »ERL :  END 


APPENDIX  C 


LANCHESTER  SIMULATION  PROGRAM  LISTING 

A  complete  listing  of  the  Lanchester  Simulation  Program  is  as  follows. 


100  'LANCHESTER  TIME  STEP  MODEL 

120  MAXFILES=2 : OPEN" LANIN" FORINPUTAS1 

121  OPEN" LANOUT"FOROUT PUT  AS2 : INPUT »1 ,NP  ,NA,ND 

122  IFNA>NDTHENMD=NAELSEMD=ND 
124  IFMD<5THENCLS: SF-1 

130  DIMAAI NA ,N0  )  ,BB1 ND ,NA  ) ,ATl NA ) ,0T( ND  > , AR( NA  )  ,DR( NO  ) 

131  DIMQA( 2 ,NA  )  ,QDt ND ) , AB<  2 ,NA  )  ,DB<  2 ,ND  )  ,SA( NA ) ,SD( ND ) ,0A(  NA ) ,0D(  ND ) 

132  'Enter  Initial  Quantities  of  Hpns,  Break  Points  And  Wpn  Types. 

134  F0RI2=1T0NA :  INPUT tfl  ,QA( 2 ,12  ) : SA< 12  )=QA(  2,12  )  :OA(  12  )=127:NEXTI2 

135  F0RI2=1T0ND :INPUT#1 >Q0( 12  )  :SD( 12 )=QDI 12  ) : 0D( 12  )  =  238:NEXTI2 

136  F0RI2=1T0NA: INPUT#1 ,ABI  1,12  ):AB(  2,12  )=AB(  1,12  )*QAi  2,12)  :NEXTI2 

137  F0RI2=1T0ND:INPUT#1,DB<  1,12): DB( 2,12  )=DB( 1,12  )*QDI  12  )  :NEXTI2 

138  F0RI2-1T0NA:  INPUT ttl  ,AT(  12  )  ‘.NEXTI2 :  F0RI2=1T0ND :  INPUTS1  ,DTl  12  )  :NEXTI2 
140  TM=0:IFSF=1THENG0SUB600 

143  F0RI1=1T0NP:PRINT#2, "STARTING  PHASE">I1 

145  'Enter  Time  Spent  In  Phase  II  and  #  of  Intervals 

146  INPUT#1,TT,NI:DT=TT/NI 

150  'Enter  Replacement  Rates  And  Attrition  Coefficient  Matrices 

152  F0RI2=1T0NA : INPUT#1 , AR< 12 ) :NEXTI2 : F0RI2=1T0ND : INPUT#1 ,DR( 12  )  :NEXTI2 

154  F0RI2=1T0NA : F0RI3=1T0ND : INPUT SI , AA(  12 ,13 ) 

156  NEXTI3:NEXTI2 

157  F0RI2=1T0ND : F0RI3=1T0NA : INPUTS1 ,BB( 12,13 ) 

159  NEXTI3:NEXTI2 

200  'Fight  Phase  II. 

202  F0RI2=1T0NI :TM=TM+DT : PRINT . 241, "Phase: ”>I1>",  Increment" >12 >"out 
of" >NI 

210  'Fight  Time  Increment  DT. 

220  'Update  number  of  attackers 

222  F0RI3  =  1T0NA:QA( 1,13  )=QA( 2,13)  :NEXTI3 : F0RI3=1T0NA : F0RI4=1T0ND 

223  ' IFDT ( 14  >  =  1THENQA( 2 ,13  )=QA( 2,13  >-AA< 13,14  )«QD( 14  )*QA( 2,13 )*DT:G0T0226 

224  'QA(  2,13 )=QA( 2,13  )-AA( 13,14  )*Q0( 14  )*DT 

225  QA(  2,13  )=QA<  2,13  )-AA(  13,14  >*<  QA(  2,13  )/QD(  14  >  >ADT(  14  >*QD<  14  >*DT 

226  NEXTI4:QA( 2,13  >=QA( 2,13 )+AR( 13  >*DT : IFSF=1THENG0SUB650 

227  NEXTI3 

230  'Update  number  of  defenders 

232  F0RI3=1T0ND:F0RI4=1T0NA 

233  'IFATI 14  )=1THENQD(  13  )=QD(  13  1-BBl  13,14  )*QD(I3  )*QA<  1,14  )*0T -.G0T0236 

234  'QD(  13  >=Q0< 13  >-BB( 13,14  )*QA< 1,14  )*DT 

235  QD(  13  )=QD<  13  )-BB(  13,14  )*(  QD(  13  )/QA(  1,14  )  JaAT(  14  )*QA(  1,14  )*DT 

236  NEXTI4:QD( 13 )=QD( 13 )+DR<  13 )*DT : IFSF=1THENG0SUB660 

237  NEXTI3 

240  G0SUB300 : NEXTI 2 

242  IFI1=NPTHENGOSUB350 : CLS: PRINT "Output  is  in  file  LANOUT .DO. " : END 
245  PRINTS2, "Status  After  Phase" >I1:G0SUB361:NEXTI1 
300  'Check  Whether  Breakpoint  is  reached. 

320  TF=0 : F0RI3  =  1T0NA : IFQAI 2 ,13  ) >AB( 2,13  )THEN325 

322  TF=1 : PRINTR2 , "Attacker  Wpn">I3>"Is  Below  Breakpoint" 

323  PRINTS2 , "  Bp  =" > : PRINTS2 ,USING"#### . ##" > AB( 2 ,13  )  > 

324  PRINT82,"  Current  Level  =" > : PRINTS2 ,USING"####.##" >QAI 2,13  ) 

325  NEXTI3 

335  FORI3=lTOND : IFQDI 13 )>0B( 2,13 UHEN340 

336  TF=1 : PRINTS2 , "Defender  Wpn">I3>"Is  Below  Breakpoint" 

337  PRINTS2,"  Bp  =" > : PRINTS2 .USING"#### . ##" >DB( 2 ,13 > > 

338  PRINTS2 , "  Current  Level  ="> :PRINT#2 .USING"#### .##" >QD< 13  ) 

340  NEXTI3 : IFTF=0THENRETURN 

350  PRINT#2,"":PRINT#2,"":PRINT#Z, "SUMMARY  AT  END  OF  BATTLE" 

351  PRINT#2,'"':PRINT#2,"Time  Elapsed  During  Battle  =”> 
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...  a  ... . 


A-V- 
Wat. 


352  PRINT#2 » USING"#### .  ##••  |  TM :  PRINT#2 , ""  :GOSUB361 
355  CLS : PRINT"0utput  is  in  file  LANOUT . DO"  (END 
361  PRINT#2,"  Att  Wpn  Breakpoint  Current  Level" 

363  F0RI3=1T0NA : PRINT#2 .USING"######" » 13 \ 

364  PRINT#! .USING"########## .##" >AB( 2,13  )  >QA( 2,13): NEXT  13 : PRINT»2 , "" 

366  PRINT#2,“  Def  Npn  Breakpoint  Current  Level" 

367  F0RI3=1T0ND :  PRINT#2  .USING"######1'  >13  \ 

368  PRINT#2 .USING"########## . ##" >DB( 2 ,13 ) jQD( 13 ) : NEXTI3 : PRINT#2 : RETURN 
600  'Set  up  output  screen 

610  PRINT"Wpn  #  Attacker  Defender" 

620  FORI 1=1T0MD: PRINTUSING"##" >11 

623  TP=2+I1*8 

625  IFI1>NATHEN630 

627  LINE( 18, TP  )-< 119,TP+4),1,B:BP=18+INT< 100*AB< 1 ,11 ) ) 

628  LINE I  BP-1, TP+1 )-( BP,TP+3  )»1,B 
630  IFI1>NDTHEN635 

632  LINE! 138,TP )-<  239.TP+4 )  ,1,B:BP=138+INT( 100*DB( 1,11 ) ) 

633  LINE! BP-1, TP+l)-< BP, TP+3),1,B 
635  NEXTI1: RETURN 

650  'Update  screen  output  of  attackers 
653  TP=3+I3*8 

655  LINE! OA( 13 ) ,TP )-( 0A(  13 )  ,TP+2  )  ,0 

656  OA( 13 )=18+INT( 100*QAl 2,13 )/SA(  13 ) ) 

657  IFOA< 13  )>118THEN0A( 13  >=118 : PRINT . ( 13440+2  ) , : G0T0659 

658  PRINT. 13*40+2,"  " 

659  LINE ( OA( 13  > ,TP  )-( 0A< 13  > ,TP+3 ) , 1 : RETURN 

660  'Update  screen  output  of  defenders 
663  TP=3+I3*8 

665  LINE ( 00 ( 13  )  ,TP  )-( 00( 13 )  ,TP+2 )  ,0 

666  0D< 13 )=138+INT( 100*QDl 13 )/S0(  13 ) > 

667  IFODl 13  »238THEN0A( 13  >=238: PRINT . 13*40+22,"*" :G0T0669 

668  PRINT. 13*40+22,"  " 

669  LINE! 0D( 13  )  ,TP  )-( 0D( 13 ) ,TP+3 ) ,1 : RETURN 


APPENDIX  D 

GEOMETRIC  PROGRAMMING  PROGRAM  LISTING 


A  complete  listing  of  the  Geometric  Programming  Program  is  as  follows. 


100  'Geometric  Programming  Program 
110  0PEN"GE0IN"F0RINPUTAS1 
120  INPUT#!, NT ,NV:K9=NT 

122  IFNT-NV<>1THENPRINT"***ERR0R:  Degree  of  Difficulty  <>  0":END 
130  INPUT#1,NC:0IMNT(NC  ) 

140  MN=0: FORI1=OTONC :INPUT#1»NT(  II  ):IFNTt  II  )>h8fTHENMN=NTI  II )  :NEXTI1 
143  DIMCTI 3 ,NC ,MN )  ,LM< NC ) ,B1<  NH-l ,NT*2 ) ,B2( NT ) ,B3(  NT ,NV ) 

145  FORI1=OTONC:FORI2=1TONTII1):INPUT#1,CT(  1,11,12 ):NEXTI2:NEXTI1 
150  F0RI1=1T0NT( 0 ) :B1(  1,11 )=1 :NEXTI1 
155  F0RI1=NT< 0 )+lT0NT:Bll 1,11 )=0 :NEXTI1 

160  F0RI1=2T0NT : F0RI2=1T0NT :INPUT#1 ,B1( 11,12 ) :B3< 12,11-1  )=B1( 11,12 ) 

162  NEXTI2:NEXTI1 

170  PRINT'"':PRINT"**COMPUTING  DELTA'S**" 

172  B2( 1  )=l:FORIl=2TONT :B2( II )=0:NEXTI1 
180  GOSUB9800 

200  CLS:I1=1: FORI2=OTONC: FOR13=1TONT( 12): CT (2, 12, 13  )=B1( 11,1 ) 

203  PRINT"DELTAl"»I2»'',"jI3»")  =  "J 

204  PRINTUSING"#####. ####"»CT< 2,12,13  ):I1=I1»1:IFI1>5THENGOSUB600 

205  NEXTI3 :NEXTI2 :GOSUB600 :CLS 

210  PRINT"" :PRINT"**COMPUTING  OPT  FN  VALUE**" 

212  FORI1=OTONC : LMl II  )=0 : F0RI2=1T0NTI II ) : LM( II  )=LMI II  )+CTI  2 ,11 ,12  ) 

214  NEXTI2:NEXTI1 
220  FS=1:FORI1=OTONC 

222  F0RI2=1T0NT< II ) : FS=FS*ICTt 1,11,12  )/CT(  2,11,12 ) )ACTI 2,11,12 ) 

224  NEXTI2:  FS=FS*(  LM(  II  )ALM<  II )  ):NEXTI1 

229  PRINT"" : PRINT"F*  s" , :PRINTUSING"####.###"»FS:GOSUB600:CLS 

230  'Compute  optimal  xln) 

232  K9=K9-1 

234  F0RI1=1T0K9: FORI 2=1T0K9: Bll 11,12 )=B3( 11,12 )  :NEXTI2:NEXTI1 

236  CC=1 : FORI1=OTONC : F0RI2=1T0NT III) 

237  CTI 3,11,12  )=C  CT(  2,11,12  )/CT(  1,11,12  )/LMIIl ) ) 

238  IFI1=0THENCT( 3,11,12  )=CT( 3,11,12  )*FS 

239  82 ( CC  )=L0G<  CTI 3 ,11 ,12  ) )  :CC=CC*1 :NEXT12 :NEXTI1 

242  PRINT"P( m, t )*  =  opt.  value  of  term  t,  constr.  m>  divided  by  its  coefficient.” 
244  FORI 1=0T0NC: FORI 2=1T0NTI II ) :PRINT"PI "  )I1}"  ,"  >I2j"  )*  =") 

246  PRINTUSING"#### . ####" iCT  I  3  , 1 1 , 1 2  ) : NEXTI 2 : GOSUB600 : NEXTI1 : CLS 
250  PRINT'"': PRINT"**  Computing  Opt  Values  Of  XI n)  **” 

260  GOSUB9800 

270  CLS:F0RI1=1T0K9:PRINT"X*I "ill j" )  =  "> 

272  PRINTUSING"###### . #####" > EXP I Bll 1 1 , 1 ) ) 

273  IFI1>5THENGOSUB600 
275  NEXTI 1:GOSU8600: END 

600  INPUT"**  Hit  ENTER  To  Continue:  "jZ9: RETURN 

9800  'Simultaneous  Linear  Equation  Subroutine:  Ax=b 

9815  'Invert  Matrix  A 

9820  F0RK7=K9+1T02*K9: F0RK8=1T0K9 

9822  IFK7=K8+K9THENB1<  K8,K7  X1ELSEB1I K8  ,K7  )=0 

9825  NEXTK8:NEXTK7 

9830  F0RK7=1T0K9 

9835  IFB1I  K7,K7  )*SGN1  Bll  K7,K7 )  X1E-8THENGOSUB9910 

9840  K2=1/B1(K7,K7  ) :  F0RK6=1T02*K9:B1(K7,K6  )=B1(  K7,K6  )*K2:NEXTK6 

9842  IFK7=K9THEN9865 

9845  F0RK8=K7+1T0K9 : IFB1I K8 ,K7  )=0THEN9860 
9850  K2=-B1(K8,K7) 

9855  F0RK6=K7T02*K9:B1IK8,K6  )=B1(K8,K6  )♦( K2*B1I  K7,K6  )  ):NEXTK6 
9860  NEXTK8:NEXTK7 
9865  F0RK7=K9T02STEP-1 
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9870  FORK8=K7-1T01STEP-1:IFB1<K8,K7)=OTHEN9885 
9875  K2=*B1<K8,K7) 

9880  F0RK6=1T02*K9: BX ( K8 ,K6 ) =B1< K8 ,K6 )♦ C K2*B1( K7  ,K6 )) : NEXTK6 

9885  NEXTK8:NEXTK7 

9890  'tfcilt  A  Inverse  by  b 

9894  F0RK7=1T0K9:B1(K7>1  )=0 : F0RK8=1T0K9:B1<  K7»l  )=B1( K7 »1  )*BHK7>K8+K9  )*B2(K8 ) 
9896  NEXTK8:NEXTK7: RETURN 
9900  ‘Error  Routine 

9903  IFERL>9700AND£RR=11THENPRINT"! ! f ERROR:  Matrix  Is  Not  Inver table ! • • " : END 
9905  PRINT“Error  Code"  jERR>‘‘In  Line" iERL:END 
9910  ‘SWITCH  ROWS 

9915  F0RK5=K7+1T0K9: IFB1< K5.K7  )»SGN<  Bl(  K5.K7 )  X1E-8THEN9940 
9920  F0RK4=1T0K9»2 : K3  =B 1 ( K7 , K4 ) : B1 1 K  7 , K4 )  =B1 ( K5 , K4 ) 

9930  B1 ( K5  »K4 ) =K3 : NEXTK4 : RETURN 

9940  NEXTK5: PRINT "Error:  Matrix  Not  Invertable" : END 


APPENDIX  E 

MATRIX  ALGEBRA  PROGRAM  LISTING 


A  complete  listing  of  the  Matrix  Algebra  Program  is  as  follows. 

100  CLS : PRINT"" : PRINT"  ***  MATRIX  ALGEBRA  PROGRAM  ***": PRINT"" 

105  PRINT"IS  INPUT  MATRIX,  'MATIN. DO'  IN  RAM?" : INPUT"  0=NO,  1=YES">FF 

107  I FFF=1THEN0PEN"MATIN"F0RINPUT AS1 

300  PRINT"**Enter  The  Single  Largest  Dimension  of" 

305  INPUT"The  Largest  Matrix  To  Be  Processed:  "jK 

310  DIMAK  3,K,K  )  »Bll K+l ,K*2  )  ,R(  4  ),C(  4 )  ,DET(  2  )  :MI  =  1 :0F=1 :SF=0 

501  CLS:EF=0:PRINT"****MATRIX  ALGEBRA  PROGRAM  MENU****" 

504  PRINT"  1.  Enter  Starting  Left  Side  Matrix" 

505  PRINT"  2.  Matrix  Inversion" 

506  PRINT"  3.  Matrix  Addition":PRINT"  4.  Matrix  Multiplication" 

508  PRINT"  5.  Simultaneous  Linear  Equations" 

509  PRINT"  6.  Print  Current  Answer  Matrix": PRINT"  7.  Other  Options"* 

510  INPUT"  **Enter  Number: "jCH 

512  IFCH=1THENMI=1:Z9=0:GOSUB7006 

513  IFCH=2THENMI=1:GOSUB2000 

514  IFCH=3THENGOSUB3000 

515  IFCH=4THENGOSUB5000 

516  IFCH=5THENGOSUB4000 

517  IFCH=6THENGOSUB6000 

518  IFCHO7THENG0T0501 

520  CLS:PRINT"***MORE  CHOICES:": PRINT"  1.  Determinant" 

524  PRINT"  2.  Matrix  Integer  Exponentiation" 

526  PRINT"  3.  Store  Current  Matrix" 

530  PRINT"  4.  Retrieve  Stored  Matrix" 

531  PRINT"  5.  Scalar  Multiplication" 

532  PRINT"  6.  Other  Options" :INPUT”**Enter  Number:  "*  CH 
540  IFCH=1THENMI=1:GOSUB1000 

548  IFCH=2THENMI=1 : GOSU87600  -  - 

549  IFCH=3THENGOSUB8000 

550  IFCH=4THENMI=1:GOSUB8200 
560  IFCH=5THENMI=1:GOSUB5100 
570  GOT0501 

700  'PAUSE  CONTROL 

702  INPUT"**  Hit  Enter  To  Continue"*Z9:RETURN 
800  'INTERMEDIATE  MODIFICATIONS 
810  PRINT''**Modify  The  2nd  Matrix?" 

812  INPUT"  0=No»  l=Invert,  2=Scalar  Multiply:  ">Z9 
815  IFZ9=0THENRETURN 

820  MI=2: IFZ9=1THENGOSUB2000E LSEGOSUB5100 
825  GOT0810 

1000  'CALC  DETERMINANT 
1005  IFR(MI)=C(MI)THEN1010 

1007  PRINT"ERROR:  Number  of  rows/columns  not  equal:” 

1008  PRINT"  MATRIX  IS  NOT  INVERTABLE  * " : GOSUB700 : EF=1 : RETURN 

1010  F0RI1=1T0R(  MI  ) :  F0RI2=1T0C(  MI )  :B1(  11,12  )=A1<  MI  ,11,12  )  :NEXTI2  -.NEXTIl 

1020  OETIMI  )=1:F0RI1=1T0R(MI ) 

1021  IFB1(  11,11  )*SGN(  Bl(  11,11 )  X1E-10THENGOSUB1900ELSE1023 

1022  IFEF=1THEN1008 

1023  DET(  MI  )=D£T<  MI  )*B1(  II ,11 ) :IFI1=RI  MI  ITHEN1080 

1025  F0RI3=1T0C(MI )  :B1( 11,13  )=B1( 11,13  )/Bl( 11,11 )  :NEXTI3 
1030  F0RI2=I1+1T0R( MI  ):IFB1( 12,11  )=0THEN1060 

1040  F0RI3=I1T0C(MI):BHI2,I3)=B1(I2,I3)-(BHI2,I1)*B1(I1,I3U:NEXTI3 
1060  NEXTI2:NEXTI1 

1080  IFOFOITHENRETURN 

1081  PRINT"**0et.  Of  Matrix  "*MI>"  Is:  "» 

1082  PRINTUSING"ltt»###. #####" *DET( MI ):GOSUB700 
1090  RETURN 

1900  'SWITCH  ROWS 


1910  F0RJ=U+1T0R<  MI ) :  IFB1U  ,11  )«SGN(  Bl(  J  ,11 )  X1E-10THEN1940 
1920  FORJlalTOCl MX )»2  :TE=B1< II ,J1 ) :Bll II, J1  )=B1( J, J1 ) 

1930  Bl(  J,J1  )=TE :NEXTJ1:GOT01950 
1940  NEXTJ:EF=1: RETURN 
1950  OET(  MI  )*-DET<  MI)  .-RETURN 
2000  'MATRIX  INVERSION 

2010  OF=0 :GOSUB1000 : IFOETI MI  )*SGN(  0ET( MI ) )>1E-100REF=1THEN2017 

2015  PRINT" TERROR:  Determinants .  MATRIX  NOT  INVERTABLE!":GOSUB700:EF=1 

2017  IFEF=1THENRETURN 

2020  F0RI1=1T0R(MI  ):F0RI2=1T0C(MI ):B1( II ,12 XAlfMI ,11,12 ):NEXTI2;NEXTI1 
2030  FORIlaCCMI  )+lT02#C( MI )  :F0RI2=1T0R( MI ) 

2032  IFI1=I2+R( MI  )THENB1( 12,11 )=1ELSEB1( 12,11 )=0 
2035  NEXTI2:NEXTI1 
2040  F0RI1=1T0C(MI  ) 

2045  IF81(  11,11  )*SGN<  Bl(  II  ,11 )  X1E-10THENGOSUB1900ELSE2055 

2046  IFEF=1THEN2015 

2055  ML=1/B1<I1,I1).-F0RI3=1T02*C<MI):B1<H,I3XB1(I1,I3)*ML:NEXTI3 
2057  IFI1=CCMI  )THEN2080 

2060  F0RI2sIl+lT0R(MI )  :IFB1(  12,11  X0THEN2075 
2065  ML=-Bll 12,11 ) 

2070  F0RI3=I1T02*C(MI  ) :B1(  12,13  >=B1<  12,13  )+CML*Bl(  11,13  ) ) :NEXTI3 

2075  NEXTI2:NEXTI1 

2080  FORI1=CCMI)T02STEP-1 

2100  F0RI2=I1-1T01STEP-1:IFB1< 12,11 X0THEN2130 
2110  ML=-B1( 12,11) 

2120  F0RI3=1T02*C(MI )  :B1(  12,13  )=B1(  12,13  )♦( ML*B1(  II  ,13  ) )  :NEXTI3 

2130  NEXTI2:NEXTI1 

2140  F0RI1=1T0C( MI  ):F0RI2=1T0R1MI  ) 

2145  All  MI  ,12,11  )=B1< I2,I1+C<  MI ) )  .-NEXTI2  .-NEXTI1 
2190  OF=l: RETURN 
3000  'MATRIX  ADDITION 

3010  MF =2 : MI s  2 : GOSUB7000 : GOSUB800 : 1 FE  F = 1THENRETURN 

3015  F0RI1=1T0R(  1 ) :  F0RI2=1T0C(  1 ) :  All  1 ,11,12  )=A1(  1,11,12  XA1I  2,11,12) 

3020  NEXTI2: NEXTI1 : GOSUB6000 :MF=0: RETURN 
4000  'SIMULTANEOUS  LINEAR  EQUATIONS 

4010  CLS:PRINT"»*Solves  Ax=b.  Choices:": PRINT"  1.  Enter  b  Vector." 

4012  PRINT"  2.  Change  An  Element  In  Matrix  A." 

4013  PRINT"  3.  Solve  Current  Ax=b." 

4014  PRINT"  4.  Return.": INPUT"  «  Select  A  Number:  ")CC 
4020  IFCC=1GOT04040 

4022  IFCC=2GOT04050 

4024  I FCC =3GOT04060  -  - 

4026  IFCC=4THEN  RETURN 
4035  GOTO  4000 

4040  MI=2:R( 2  )=C( 1  ):C<  2  )=1:GOSUB7040:GOT04000 

4050  INPUT "**Row ,  Column  Of  Matrix  A  To  Be  Changed:  "jRD.CO 

4052  PRINT"  -  Enter  Row"iRD|",  Column" »C0 )":'*)  :INPUTA1(  1,RD,CD  ) :GOT04000 

4060  MI=1:SF=1;OF=0.-GOSUB8000:GOSUB2000 

4064  IFEF=0THEN4070 

4065  PRINT"*Solution  Not  Uniquely  Determinable" :GOSUB700: RETURN 
4070  GOSUB5000:CC=0:FORIl=lTOR(2):CC=CC+l:PRINT"x("vIl»")  =  ") 

4075  PRINTUSING"##### . ####"  )B1(  II ,  1 ) :  IFCO6THENG0SUB700 : CCS 
4080  NEXTI1 : GOSUB700 :SFS: G0SUB8200 .- GOT04000 

5000  'MATRIX  MULT 

5010  MF=1:IF  SF=1THEN5020 

5015  MIa2:GOSUB7000:GOSUB800 : IFEF=1THENRETURN 

5020  R(  4  XR(  1  ):C(  4  )=C(  2  ) : F0RI1=1T0R( 4 ): F0RI2=1T0C( 4 )  :B1(  11,12  )S 
5022  F0RI3*1T0C< 1 ) :B1( 11,12 )SA1( 1,11,13  )»A1( 2,13,12  )+Bl( 11,12  ) 

5024  NEXTI3:NEXTI2:NEXTI1:MFS 
5050  IF  SFSTHENGOSUB7500.-GOSUB6000 
5060  RETURN 
5100  'SCALER  MULT 

5110  INPUT-'Enter  Scalar  Multiplier :" )SM: F0RI1=1T0R( MI ): F0RI2=1T0CI MI ) 
5115  A1IMI,I1,I2XA1(MI,I1,I2  )*SM :NEXTI 2 :NEXTI1 : RETURN 
6000  'PRINT  OUTPUT  MATRIX 

6010  PRINT"  **  Current  Answer  Matrix:":CCS:  F0RI1=1T0R<  1  ):CCSC+1 
6012  F0RI2*1T0C< 1 ): PRINTUSING"####.####" )A1(1,I1,I2)) :NEXTI2:PRINT"" 

6050  IFCC*5THENGOSUB700:CCS 
6060  NEXT1 1 ; GOSUB700 : RETURN 


7000  'MATRIX  INPUT 

7001  CLS : PRINT"" : PRINT "Will  This  Matrix  Bo  On:" 

7002  INPUT"  0=Left,  l=Right"jZ9:  IFZ9=1THEN7006 

7003  R( 2  )=R( 1 ) :C( 2  )*C( 1 ) : F0RI1=1T0R< 1 ) : F0RI2=1T0C( 1 ) 

7004  Al( 2,11,12  >=A1< 1,11,12  ):NEXTI2:NEXTI1:MI=1 

7006  CLS:IFMI=2THEN7008 

7007  PRINT"**Choices  For  Loft  Hand  Matrix" :GOT07009 

7008  PRINT "**Choices  For  Right  Hand  Matrix:" 

7009  PRINT"  1.  Enter  Matrix  From  Keyboard" 

7010  PRINT"  2.  Retrieve  Stored  Matrix" :IFFF<>1THEN7012 

7011  IFE0FI 1  )=0THENPRINT"  3.  Enter  Matrix  From  MATIN. DO" 

7012  INPUT"#*Enter  A  Number:  " iZ9:IFZ9=lTHEN7015 

7013  IFZ9=3THEN7018 

7014  R(  MI )=R( 3 ) :C( MI J=C« 3 )  :GOT07020 

7015  PRINT"  *»Enter  The  Rows,  Columns" 

7017  INPUT "In  The  Next  Matrix:  "  JRIMI )  ,C<MI  ):GOT07020 

7018  INPUT#1,R(MI ),C(MI ) 

7020  IF  MFO1THEN7030 

7021  IFR( 2  )=C< 1 ITHEN7030 

7022  PRINT“**ERROR:  Columns  in  LEFT  MATRIX  ="iCU) 

7024  PRINT"  Rows  In  Right  Matrix  =">R(2) 

7026  PRINT'These  Must  Be  Equal  For  Matrix  Mult !!" :G0SUB 700 : EF=1 :GOT07006 

7030  IF  MFO2THEN7035 

7031  IF(  R(  1  )=R(  2  IANDCI 1 1  =C C  2  )  1THEN  7035 

7032  PRINT“**ERR0R : D imens ions  For  Both  Input" 

7034  PRINT "Mat rices  Must  Be  Equal! !":GOSUB700:EF=1:GOT07006 

7035  IFZ9=2THENGOSUB8200: RETURN 

7036  I FZ9=3THEN7050 

7037  PRINT"  **Fill  Matrix  Row  By  Row: ":PRINT"" 

7040  F0RI1-1T0RI MI  ):F0RI2=1T0C(MI ) 

7042  PRINT"-Enter  Row">Il)"And  Column"}I2j":") 

7044  INPUTAlt MI ,11, 12 ):NEXTI2: PRINT"" :NEXTI1: RETURN 
7050  F0RI1=1T0R(MI ) : F0RI2=1T0C< MI  ) :  INPUTR1  ,A1(  MI  ,11 ,12  ) 

7052  NEXT  T  2 : NEXTI 1 : RETURN 

7500  'MAKE  ANSWER  MATRIX  THE  FIRST  MATRIX  FOR  THE  NEXT  OPERATION 
7510  R(  1 )=R( 4 ):C( 1 )=C(4 1: F0RI1=1T0RC 1 ) : FURI2=1T0C( 1 ) 

7512  Al( 1,11,12  )=B1( 11,12 ): NEXTI 2: NEXTI 1: RETURN 
7600  'MATRIX  INTEGER  EXPONENTIATION 

7610  CLS:PRlNT"":lNPUT"*»Enter  Positive  Integer  Exponent:  "jXP 
7620  R<  2  )=R( 1 )  :C( 2  l=C( 1 1 : F0RI1=1T0R( 1 1 : F0RI2=1T0C(  2 ) 

7622  All  2,11,12 )=A1( 1,11,12 ):NEXTI2:NEXTI1 

7630  SF=1: F0REX=2T0XP : GOSU85020 : GOSUB7500 : NEXTEX : GOSUB6000 : SF=0 : RETURN 
8000  'STORE  AK1,,) 

8010  R<  3  )=R( 1 1  :C( 3  )=C< 1 ) : F0RI1=1T0R<  3  ) : F0RI2=1T0C( 3  ) 

8012  All  3,11,12  )=A1( 1, II, 12) : NEXTI 2: NEXTI 1: RETURN 
8200  'RETRIEVE  THE  STORED  MATRIX 

8210  R( MI  )=R«  3  )  :C<  MI  )=CI 3  ) : F0RI1=1T0R( MI ) : F0RI2=1T0C( MI ) 

8212  A1(MI,I1,I2 )=All 3,11,12 ): NEXTI 2: NEXTI 1 : RETURN 


APPENDIX  F 

NUMERICAL  INTEGRATION  PROGRAM  LISTING 


A  complete  listing  of  the  Numerical  Integration  Program  is  as  follows. 


1200  'Numberical  Double  Integration: Steven  H.  Cary: 24  Apr  86 

1201  DIMA2( 6 >6  ):TL=.001 

1205  CLS : PRINT’"' : PRINT"  **  Double  Integration  **" 

1206  PRINT"  Romberg  Algorithm" 

1210  PRINT"0=Edit  Fisiction  To  Be  Integrated." 

1211  PRINT"l=Edit  Limits  Of  Integration. " 

1213  PRINT"2=Edit  Tolerance}  Current  Tol.s"}TL 

1215  PRINT"3=Calculate  Integral.  ":INPUT”£nter  0,  1,  2,  or  3:*'»Z9 

1216  IFZ9=0THENEDIT1285-1288 

1217  IFZ9=1THENEDIT1291-1298 

1218  IFZ9=2THENINPUT"Tolerance  =  " jTL:GOTQ1205 

1220  CLS : PRINT"" : PRINT"  !!Be  Patient! !":PRINT"" 

1230  N=2:G0SUB1293:DX=(XU-XL)/2 

1240  F0RJ9=1T06:DX=DX/2:N=N*2 

1242  X=XU:GOSUB1296:GOSUB1280:A2( J9»l )=SS*DY 

1245  X=XL:G0SUB1296 : GOSUB1280 : A2< J9 , 1  )=A2t  J9 ,1  )»SS*0Y 

1250  FORJ8=2TON:X=X+DX:GOSUB1296:GOSUB1280 

1251  A2(J9,11=A2< J9,1>+2»SS*0Y:NEXTJ8 

1252  A2U9,l)=A2(J9,l)*0X/3 
1255  IFJ9=1THENNEXTJ9 

1260  F0RJ8=1T0J9-1 

1261  AZl  J9, J8*l  )=A2I  J9, J8)+l  J  A2I J9,  J8  )-A2l  J9-1, J8 ) )/( 4A  J8-1 ) )  :NEXTJ8 

1262  Tl=A2t J9iJ9  !-A2<  J9 , J9-1 ) : IFSGNI T1  )*T1-TL>0THENNEXTJ9ELSE1264 

1263  PRINT"Tolerance  of"»TL»"not  met  after  five  extrapolations" 

1264  IN=A2(J9,J9) 

1265  PRINT "Integral  =•'} : PRINTUSING"###### . #####”» IN 

1266  PRINT"  Actual  Tolerance  =  "> : PRINTUSING"*# . ####"»T1*SGN(T1 ) 

1267  SOUND 1567 , 10 : S0UND1244 , 10 : SOUND1046 , 10 : SOUND 783 , 20 

1268  SOUND1046  »10 :SOUND783  >40 

1269  INPUT"Hit  Enter  To  Continue: "}Z9:GOT0120S 

1275  F0RJ7=1T06 : F0RJ6=1T0J7 : PRINTUSING"## . ###" > A2I J7 , J6  ) } 

1276  NEXTJ6 : PRINT"" : NEXTJ7: INPUTZ9 : RETURN 

1280  REM  Simpson's  Rule  Sum 

1281  Y=YU : GOSUB1285 : SS=F : Y=YL : G0SUB1 285 : CS=SS+F 

1282  FORJ5=2TOl N/2 ) :Y=Y+DY : G0SUB1285 : SS=SS+4*F : Y=Y«OY : GOSUB1285 

1283  SS=SS+2*F :NEXTJ5 : Y=Y+DY : G0SUB1285 : SS=SS+4«F : RETURN 

1285  'F(x>y)  to  be  integrated: 

1286  F=1 

1288  'X  4  Y  *  independent  variables.  Hit  F8,  then  F4  ttien  Done. 

1289  RETURN 

1290  'Limits  of  Integration: 

1291  'XLOWER/XUPPER  are  constants. 

1292  'YUPPER  4  Y LOWER  may  be  constants  or  given  in  terms  of  X. 

1293  XUPPER=1. 5707963 

1294  XL0WER=0 

1295  RETURN 

1296  YUPPER=SIN(X) 

1297  YLOHER=0 

1298  ‘Hit  F8>  Then  F4  When  Done 

1299  DY  =  (YU-YL)/(N+1):  RETURN 
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